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ABSTRACT 

The estimated stellar masses of galaxies are widely used to characterize how the galaxy pop- 
ulation evolves over cosmic time. If stellar masses can be estimated in a robust manner, free 
from any bias, global diagnostics such as the stellar mass function can be used to constrain 
the physics of galaxy formation. We explore how galaxy stellar masses, estimated by fitting 
broad-band spectral energy distributions (SEDs) with stellar population models, can be biased 
as a result of commonly adopted assumptions for the star-formation and chemical enrichment 
histories, recycled fractions and dust attenuation curves of galaxies. We apply the observa- 
tional technique of broad-band SED fitting to model galaxy SEDs calculated by the theoreti- 
cal galaxy formation model GALFORM, isolating the effect of each of these assumptions. We 
find that, averaged over the entire galaxy population, the common assumption of exponen- 
tially declining star-formation histories does not adversely affect stellar mass estimation. We 
show that fixing the metallicity in SED fitting or using sparsely sampled metallicity grids can 
introduce mass dependent systematics into stellar mass estimates. We find that the common 
assumption of a star-dust geometry corresponding to a uniform foreground dust screen can 
cause the stellar masses of dusty model galaxies to be significantly underestimated. Finally, 
we show that stellar mass functions recovered by applying SED fitting to model galaxies at 
high redshift can differ significantly in both shape and normalization from the intrinsic mass 
functions predicted by a given model. Given these differences, our methodology of using stel- 
lar masses estimated from model galaxy SEDs offers a new, self-consistent way to compare 
model predictions with observations. We conclude that great care should be taken when com- 
paring theoretical galaxy formation models to observational results based on the estimated 
stellar masses of high redshift galaxies. 
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1 INTRODUCTION 

A successful theory of galaxy formation is essential for accurately 
connecting any underlying cosmological framework with the ob- 
servable Universe. The vast dynamic range and overall complex- 
ity involved in the interplay between gas, stars and dark matter 
in galaxies strongly restricts what any model of galaxy formation 
can predict a priori. To make progress, it is necessary to use ob- 
servational results to constrain galaxy formation models. Estimat- 
ing the stellar masses of galaxies offers, in principle, a powerful 
method to characterize the galaxy population that can be compared 
directly to theoretical predictions. Unlike directly observable quan- 
tities, stellar mass is a derived quantity that can only be estimated 
from observational data through the application of a series of mod- 
els and assumptions. It is therefore critical to understand how these 
assumptions affect the reliability of stellar mass estimates and how 
any uncertainties affect global diagnostics of the galaxy population, 
such as the stellar mass function. 

* E-mail: peter.mitchell@durham.ac.uk 



The traditional approach for constraining parameters in theo- 
retical galaxy formation models is to use directly observable prop- 
erties of the galaxy population such as the luminosity function 
or the Tully-Fisher relation. This requires that the intrinsic phys- 
ical properties of model galaxies calculated using either hydrody- 
namical simulations or semi-analytic models (SAMs) can be con- 
verted into directly observable quantities. Stellar population syn- 
thesis (SPS) models (e.g. |Bruzual & Charlot|2003||Maraston|2005| > 
are combined with predicted star-formation and chemical enrich- 
ment histories of model galaxies to produce spectral energy distri- 
butions (SEDs) that can be compared with observational data. Un- 
certainties in the form of the initial mass function (IMF) of stars, 
the accuracy of SPS models and the manner in which dust attenu- 
ates the light emitted by stars make this process challenging (e.g. 
|Conroy eTa l, 2009 2010b). It can be difficult in some cases to be 
confident whether the comparison between model predictions and 
observational results does in fact show if a given model is success- 
ful in describing the underlying physics. 

The estimation of the stellar masses of galaxies from obser- 
vational data is typically achieved using the technique of broad- 
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band SED fitting, which inverts the process of generating observ- 
ables from intrinsic galaxy properties. The popularity of this tech- 
nique can be attributed to the success of multi-wavelength surveys 
such as the Sloan Digital Sky Survey (SDSS; |York et al.||2000]>, 
the Great Observatories Origins Deep Survey (GOODS; |Giavalisco| 
et al.|| 2004l> and the Cosmological Evolution Survey (COSMOS; 



Scovi 



le et al.|2007) . These surveys quantify how the galaxy pop 



ulation evolves with look-back time by utilizing broad-band pho- 
tometry from the UV to the NIR to obtain accurate photometric 
redshifts for much larger galaxy samples than would be possible 
using spectroscopy. The information and methods used to obtain 
accurate photometric redshifts by fitting stellar population models 
can readily be extended to also estimate stellar masses and other 
galaxy properties. Consequently, it has become standard practice 
to estimate these quantities whenever multi-wavelength photome- 
try is available. 

The same uncertainties associated with converting the intrinsic 
properties of model galaxies into observables also affect the accu- 
racy of SED fitting applied to observational data. However, there 
are a number of additional assumptions that must be made to es- 
timate the stellar masses of galaxies from observations. For exam- 
ple, unlike for the case of model galaxies, the star formation histo- 
ries (SFHs) of real galaxies are not known. Instead, a prior for the 
SFHs of galaxies must be adopted. Various studies have attempted 
to explore how stellar mass estimates depend on the different as- 
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eral consensus that stellar mass estimates are more reliable com- 
pared to other quantities such as star-formation rates (SFRs) and 
galaxy ages that can be estimated using the same process. How- 
ever, the reported level of uncertainty on stellar mass estimates can 
vary strongly, depending on the specific galaxy samples consid- 
ered. |CoMoyetaL]((2009} find that the stellar mass-to-light (M/L) 
ratios of star-forming galaxies at z = can only be constrained 
to within 0.3 dex at a 95% level of confidence when various un- 
certainties in SPS modelling are taken into account. This does not 
include the uncertainty associated with the choice of the IMF. Gal- 
|lazzi & Bell (2009) find that, ignoring the uncertainties associated 
with SPS modelling and dust attenuation, it is possible to constrain 
M/L ratios of galaxies with smooth SFHs to within 0.1 dex us- 
ing spectral features or a single optical colour (see also [Wilkins 
|et al.|20L3] >. |Longhetti & Sa racco (2009) consider the accuracy of 
stellar mass estimates of early-type galaxies, finding that the true 
stellar masses of mock galaxies can only be recovered to within 
a factor of « 0.3 — 0.5 dex, given the variations between differ- 
ent SPS models and metallicities. Marchesini et al. (2009) quantify 
how uncertainties associated with the assumptions made in SED fit- 
ting contribute to the error budget of the stellar mass function. They 
find that potential systematic errors associated with these assump- 
tions can dominate over other error sources, such as photometric 
redshifts or galaxy photometry. They also find that the shape of the 
mass function, particularly at the low-mass end, is sensitive to, for 
example, the assumed metallicity and the adopted dust attenuation 
law. 

We explore this topic from an alternative angle by applying 
the methods used in observational studies to estimate stellar masses 
from SEDs output by the semi-analytic galaxy formation model 
GALFORM ( |Cole et al.|2000} . We focus on understanding the vari- 
ous random and systematic errors encountered when estimating the 
stellar masses of individual galaxies and also study how these trans- 



late into errors in the stellar mass function. This exercise serves as 
a useful example of how the process of converting between observ- 
ables and intrinsic galaxy properties can have an impact on global 
diagnostics of the galaxy population. Several studies have adopted 
a similar approach by combining SED fitting with galaxy formation 
models, aiming to understand the accuracy of SED fitting in differ- 
ent scenarios <|Lee et al.|2009||Wuyts et al.|2009||Lee et aL|20T0| 
|Pforr et al.|2012| >. This method is useful because outside of the lim- 
ited number of cases where additional data is available, it is difficult 
to test the accuracy of quantities estimated using SED fitting. Fit- 
ting mock galaxy SEDs provides a means to do this and theoretical 
galaxy formation models are, in principle, a useful tool for pro- 
ducing samples of mock galaxies with realistic star-formation and 
chemical enrichment histories. In the case of SAMs, model galaxy 
samples that represent the entire galaxy population over a range of 
redshifts can be generated rapidly. This allows us to isolate and un- 
derstand different effects by considering variants of the underlying 
model. 

Lee et al. j2009 1 use model galaxy SEDs from the |Somerville| 



|& Primack| \\999) SAM to explore the accuracy of SED fitting 
in recovering the physical properties of Lyman-break galaxies 
(LBGs). They find that stellar masses are, on average, well re- 
covered for LBGs as they are underestimated, on average, by less 
than 0.1 dex. They attribute this success to two competing factors. 
Younger stars can mask the presence of older stars in LBG SEDs 
in some cases, leading to an underestimate of the total stellar mass. 
However, they also find that there is a tendency for SED fitting to 
overestimate the age of LBGs. This typically leads to overestimat- 
ing the stellar masses of some galaxies. 



Pforr et al.|f2012 1 combine SED fitting with the GallCS SAM 



(Hatton et al. 2003). They find that, in general, stellar masses 
are slightly underestimated when using standard exponentially de- 
clining SFHs in the SED fitting process. They demonstrate that 
this problem can be resolved by adopting exponentially increasing 
SFHs for star-forming galaxies. They conclude that stellar masses 
are recovered almost perfectly at redshifts z £ 2,3, where the al- 
lowed distribution of galaxy ages is fairly narrow. At lower red- 
shifts, however, they find the stellar masses of star-forming galaxies 
can underestimated by up to 0.6 dex as a result of discrepancies be- 
tween the estimated and true SFHs. They explain that this is caused 
by the larger range of possible ages at low redshift combined with 
degeneracies between age and dust. Finally, they also show that this 
problem can be circumvented by choosing to ignore dust redden- 
ing when estimating the stellar masses of star-forming galaxies at 
low redshift. This prevents the SED fitting procedure from fitting 
an unrealistically young age coupled with a large amount of dust 
reddening. 

A key difference in our methodology compared to that o f|Pforr| 
|et al.| and |Lee et ahj is that we use a physically motivated model for 
attenuation by dust. |Pforr et al.| and |Lee et al.| instead adopt empir- 
ical dust attenuation laws to calculate model galaxy SEDs, corre- 
sponding physically to a uniform foreground screen of dust placed 
between the galaxy and an observer. The difference between a fore- 
ground dust screen model and the physically motivated radiative 
transfer calculation performed in our analysis turns out to be very 
significant for our results on stellar mass estimation. It should also 
be noted that it is not our intent to follow Pforr et al. (2 012| > in 
attempting to quantify the exact level of random and systematic er- 
rors on stellar mass estimates for an exhaustive range of possible 
SED fitting configurations and filter sets. This is because any quan- 
titative results derived from the approach of combining SED fitting 
with theoretical models may be sensitive to the degree to which 
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a given model can represent the true galaxy population. Instead, 
we attempt to provide a detailed explanation of the different error 
sources we encounter when considering the overall galaxy popula- 
tion over a wide range of redshifts. This is achieved by isolating the 
different factors responsible for biasing stellar mass estimates in a 
specific set of idealized test cases. 

The outline of this paper is as follows. In Section [2] we in- 
troduce broad-band SED fitting, explain some of the underlying 
assumptions that are involved in the process and outline the param- 
eter choices we make in this study. Section|3]gives a brief overview 
of GALFORM and the specific models which we use in this study. 
We also explain how intrinsic galaxy properties are transformed 
into observables in the context of the assumptions made in SED 
fitting. In Section [4] we present results of performing SED fitting 
on model galaxies, focusing on exploring the systematics that affect 
the recovery of the stellar masses of individual galaxies. We present 
results for the stellar mass function over a range of redshifts in Sec- 
tion]^] Finally, we discuss and summarize our results in Section|6] 
and Section|7] 



2 BROAD-BAND SED FITTING 

We seek to understand the relationship between the stellar mass es- 
timated from observations using SED fitting and the stellar mass 
predicted by theoretical models. Instead of using broad-band pho- 
tometry from observed galaxies, we fit the broad-band magnitudes 
of model galaxies predicted by the semi-analytic model GALFORM. 
The precise details of the method used to perform SED fitting in 
different observational studies typically vary very little. Detailed 

BoFl 



descriptions and discussion of the method can be found in 



zonella et al. 



Taylor et al. 



l [2000l ), [S^mTetli^ ( [2007l >, |Walcher et aTJpOTT 



and 

201 \) . In this section we provide an overview of SED 
fitting as a method of estimating stellar mass. We also describe our 
parameter choices for the SED fitting procedure used in this study. 



2.1 Overview 

Broad-band SED fitting works by comparing a grid of tem- 
plate galaxy SEDs to observational data. Typically, a maximum- 
likelihood method is then used to decide which template best fits 
the data (although see |Taylor et al.|2011| |Salim et al.|2007| for a 
discussion of alternative statistical techniques). This is achieved by 
first minimizing x f° r eac h template SED, then choosing the tem- 
plate SED with the smallest associated \ 2 value. This corresponds 
to choosing the mode of the likelihood distribution. \ 2 is calculated 
by summing over all available photometric bands bands using 



x 2 = E 



-^galaxy,™ S -ftemp,n 
17 n 



(1) 



where Fgaiaxy.n and F tenl p,n are the fluxes of the galaxy and tem- 
plate, respectively, in the nth band, s is a normalization factor and 
er n is the la flux error associated with a galaxy in the nth band. 
The normalization factor s is calculated such that \ 2 is minimized 
for each template SED using 
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where a choice can be made regarding which bands are included in 
the summation. We choose to follow standard practice by simply 
summing over all available photometric bands, as in Eq. [T] This 
process yields a best-fitting M/L ratio for each galaxy which can 
then be used to estimate the stellar mass of a given galaxy. 

Template galaxies SEDs are generated using publicly avail- 
able SPS models (e.g. |Bruzual & Charlot||2003| |Maraston|[2005l 
|Conroy et al.|2009) . SPS models predict the spectra of simple stel- 
lar populations (SSPs), a group of stars with the same age and 
metallicity and a distribution of initial masses given by the stellar 
initial mass function (IMF). These SSP spectra are then convolved 
with an assumed parametrization for the typical SFH of a galaxy. It 
is well established that various degeneracies make it very difficult 
to place strong constraints on galaxy SFHs from photometric data 
alone, unless strong priors are adopted (e.g. Ma raston et al.|20l0| >. 
It is therefore standard practice to assume a simple parametrization 
for the SFH which can represent a broad range of possible SFHs 
without creating a prohibitively large parameter space over which 
to search. By far the most common choice of parametrization used 
for the SFH of low and intermediate redshift galaxies is an expo- 
nentially declining SFH. It should be noted, however, that there 
are numerous studies which have advocated alternatives, particu- 
larly for high redshift galaxies (e.g. |Lee et a l. 2010, Micha lowski| 
|et al.|2 012; Pacifici et al. 2013). The exponentially declining SFH 
is parametrized by the time since the onset of star-formation, t age , 
and the e-folding timescale, r. SPS models also output the mass 
returned from a SSP back into the interstellar medium (ISM) as 
a function of age, which in turn is used to predict the M/L ratio 
of each template galaxy. To reduce the size of the parameter space 
that needs to be searched over, it is typically assumed that all of the 
stars in each template galaxy have a single stellar metallicity, Z*. 
Furthermore, the number of metallicity points in the parameter grid 
is usually very small due to the sparse metallicity grid made avail- 
able for publicly available SPS models. For example, there are only 
5 metallicities available for the |Bruzual & Charlot| ( |2003l (BC03) 
SPS model. Unlike what is done in theoretical models, it is not stan- 
dard practice to interpolate between metallicities in SED fitting. 



2.2 Dust attenuation 

Observed galaxy SEDs are a product of the intrinsic galaxy SED 
produced by stellar emission which is then attenuated according 
to radiative transfer through intervening gas and dust. It is stan- 
dard practice in SED fitting to account for absorption by neutral 
hydrogen in the IGM using the prescription from M adau|fl995| >. 
This is adopted in both GALFORM and in the SED fitting proce- 
dure used in this study. Attenuation by dust in the ISM of galaxies 
is a substantially more complex radiative transfer problem. Dust in 
galaxies can be concentrated around young stars or distributed dif- 
fusely throughout regions of the ISM. Given the lack of information 
available on the relative spatial distribution of stars and dust in dis- 
tant galaxies, attenuation by dust it usually accounted for in SED 
fitting using the empirical Calzetti dust attenuation law (Calzetti 
|et al.|1 994 2000). The Calzetti attenuation law has a fixed shape 
which is different from the dust extinction curve in the local ISM. 
The star-dust geometry that is implicitly assumed when applying 
the Calzetti law corresponds physically to a uniform, foreground 
dust screen placed between the observer and the stellar populations 
in a given galaxy. Making this assumption has the advantage for 
SED fitting in that the Calzetti is consequently only characterized 
by only a single parameter, the reddening E(B — V), defined as the 
difference between the observed and intrinsic B — V colour. The 
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shape of the Calzetti law was derived from a sample of 30 local 
starbursts (Calzetti et al. 1994) and the normalization was derived 
from a sub-sample of only 4 local starbursts ( Calzetti et al. 2000). 

A number of studies have attempted to assess how well the 
Calzetti law can reproduce the attenuation curves of different 
galaxy types. At low redshift, Wild et al.| ( |2011| > apply a pair- 
matching technique to study the shape of the attenuation curves 
of SDSS spirals. They find evidence that the shape of the optical 
dust attenuation curves of local star-forming galaxies is strongly 
dependent on galaxy inclination. Specifically, they show that face- 
on spirals have steeper optical attenuation curves than the Calzetti 
law, whereas edge-on spirals are consistent with the Calzetti law. 
They also find that the slopes of the near-infrared (NIR) attenuation 
curves are consistent with Milky-Way extinction or Calzetti law at- 
tenuation curves, independent of inclination. In the UV, they find 
evidence for a bump in the attenuation curve of spirals at 2175 A 
(see also Conroy et al. 2010a). This feature is absent from the 
Calzetti law and could have a significant impact on the interpre- 
tation of the properties of high redshift galaxies (Gonzal ez-Perez| 
|etal.|2013| >. 

Another method that is used to investigate the dust attenua- 
tion properties of galaxies involves measuring how the ratio of far- 
infrared (FIR) to ultraviolet (UV) flux, IRX, varies as a function of 
the UV spectral slope, j3 lBell|2002||Goldader et al.|2002[|HoweIIl 



|et al.|2010| [Murphy et al.|20li||Penner et al.|2012| e.g.). The po- 
sition of galaxies on the IRX-/? plane can then be compared with 
the relationship derived for local starbursts JMeurer et al.||1999| 
|Overzier et al.|20fl) . The |Meurer et al.| < |1999| > relation was derived 
from the same galaxy sample used to derive the Calzetti law and so 
this comparison tests whether the Calzetti law is applicable objects 
other than modestly starbursting local galaxies. Bell (2002) were 
the first to show that local star-forming galaxies lie below this re- 
lation such that there is less UV attenuation for a given value of f3. 
Goldader et al. ( 2002 ) were the first to show that local ultralumi- 
nous infrared galaxies lie above this relation such that the Calzetti 
law underestimates the total UV attenuation for these objects. 

Various studies have also applied this method to high redshift 
star- forming galaxy sam ples (e.g. |Murphy et al.|2011| |Buat et aT] 
[20T2{[Reddy et al.|2012]|Penner et al.|2012) . This exercise is dif- 
ficult because FIR SEDs are typically available only for the most 
extreme dusty galaxies at higher redshifts, although a stacking anal- 
ysis can ameliorate this problem (Reddy et al. 2012). There is ev- 
idence for consistency between the high redshift and local IRX-/3 
relations (e.g. Reddy et al. 2010 2012 1, although it has also been ar- 
gued that specific object classes can be offset from the local relation 
(e.g. |Murphy et al.|201 l||Penner et al.|2012| >. |Buat et al.| ( [20T2) ap- 
ply a different approach and analyse a sample of 75 1 UV-selected 
galaxies at z € 1,2, fitting the full UV-FIR SEDs using a SED 
fitting procedure that features a generalized form of the Calzetti 
law ( |Noll et al.|20 09 1. They find evidence for a steeper attenuation 
curve in the UV than the canonical Calzetti law for 20% of their 
sample, and also a UV bump. 

2.3 Filter and parameter choices 

The final step in producing template SEDs involves convolving 
with the broad-band filters used in given observational dataset. 
Deep multi-wavelength surveys such as GOODS have many photo- 
metric bands available, spanning all the way from the UV through 
to the radio. Wide-area surveys on the other hand such as SDSS 
typically have only optical broad-band photometry available. For 
simplicity, we use a single filter set across a wide range of redshifts 



with the exception of Section [5~T| where we consider LBG samples. 
Filters blueward of the Lyman limit at 912 A in the rest frame are 
excluded from the fitting process. We do not include artificial red- 
shift or flux errors, setting a n in Eq.[T]to 10% of the model galaxy 
flux ^galaxy,™ for each band. These error sources become impor- 
tant at high redshift but can be understood without the need for a 
theoretical model and are not particularly relevant for understand- 
ing errors in stellar mass estimates associated with the assumptions 
made in SED fitting. However, ignoring them entirely means that 
any quantification of the errors in stellar mass estimates given in 
this paper should be considered as lower limits. We do consider the 
effect of artificially perturbing model fluxes when exploring LBG 
samples in Section [BTT] where it becomes important to include de- 
tection criteria in order to robustly compare model predictions with 
observational data. 

In this study we use two filter and SED fitting parameter sets, 
as outlined in TableQ] A common feature of both sets is that we use 
|Bruzual & CharlotlpOOl] ) SPS models and the Calzetti law for SED 
fitting. t agc is always constrained to be less than the age of the Uni- 
verse at the given redshift. We also follow Santini et al. (2012]) and 
others by excluding super-solar metallicity templates for z ^ 1. 
Note that this is not enforced for model galaxies in GALFORM. 
We refer to the first parameter set as the standard parameter grid 
because it is designed to be broadly representative of the choices 
made in observational studies of low to intermediate redshift galax- 
ies where Spitzer IRAC imaging between 3.6/^m and 8fim is often 
available (e.g. |Ilbert et al.|2010||Santini et al.|2012||Morflock et al.| 
|201 1\ . It uses the exponentially declining SFH typically used for 
galaxies at low and intermediate redshift and is therefore charac- 
terized by t ago , r, Z t and E(B — V). We use a Salpeter IMF for 
this parameter set despite the fact that the GALFORM models we 
consider typically use a Kennicutt IMF to demonstrate how the sys- 
tematic uncertainty on the IMF compares against other sources of 
error in stellar mass estimation. We modify this choice of IMF in 
the templates to a Chabrier IMF in Section [5] where we consider 
model predictions of the stellar mass function. This choice is made 
so that the stellar masses of model galaxies estimated from SED 
fitting are consistent with the observational studies with which we 
compare. 

We also use a second parameter set, deliberately constructed 
to closely resemble the choices made by |Lee et al.| ( |2012") , who use 
SED fitting to estimate the stellar masses of LBGs at z = 4 and 
z — 5. We refer to this parameter set as the LBG parameter grid 
and use it in Section [5~T| It uses an exponentially declining SFH, 
including the limit of r — > oo, corresponding to a constant SFH. 
When considering LBGs, it is important to consider both the effect 
of photometric errors and non-detections where galaxies drop be- 
low the sensitivity limit of the survey. Instead of setting a n in Eq.[T] 
to 10% of the model galaxy flux -Fgaiaxy.n, in this case we use the 
5-sigma limiting magnitudes listed in Table 1 of |Lee et aL|{2012| > to 
determine a n . We also consider the effect of artificially perturbing 
fluxes using a Gaussian distribution with a again taken from Table 
1 in |Lee et al.|{2012| l. To facilitate a self-consistent comparison, 
we follow the same procedure for non-detections described in Lee 
|et al.| i 2012[ ), where non-detected bands are used as upper limits 
in Eq. T| if s Jt e m P ,n exceeds the upper limit in that band. We use 
the LBG dropout selection criteria (both colour and S/N selection) 
given by Equations 1-13 in |Stark et al.| ( |2009) . These criteria select 
B435, Veo6 and 4775 dropouts to create samples of LBGs at z — 4, 5 
and 6 respectively. We use a Chabrier IMF for this parameter set. 
The allowed parameter values for both parameter sets are listed in 
Table [I] 
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Standard Parameter Grid 



Filters B435, V 60 6, R, «775, 2850, J, H, K, 3.6, 4.5, 5.8, 8.0/Ltm 

IMF Salpeter 

tagc/Gyr 0.1, 0.11, 0.13, 0.14, 0.16, 0.18, 0.2, 0.23, ... 

r/Gyr 0.1, 0.3 0.6, 1, 2, 3, 4, 5 7, 9, 13, 15, 30 

Z*/Z 2.5,1,0.4,0.2,0.02 

£(B - V) 0, 0.03, 0.06, 0.1, 0.15, 0.2, 0.25, 0.3, ... 1 

Lyman-Break Galaxy Parameter Grid 



Filters 
IMF 

*age/Gyr 
r/Gyr 
Z*/Z 



B435. V606. «775. -2850. J, H, K, 3.6, 4.5, 5.8/im 
Chabrier 

0.1, 0.11, 0.13, 0.14, 0.16, 0.18, 0.2, 0.23, ... 

0.1,0.2 0.3,0.4, 0.6, 0.8, 1.0, 00 

1,0.2 



E(B - V) 0, 0.025, 0.05, 0.075 ... 0.95 

Table 1. Parameter grids for SED fitting. The top section outlines the stan- 
dard parameter grid we use for the majority of our results. The bottom sec- 
tion outlines the parameter grid used in Section [5~T] for exploring LBG sam- 
ple selection. i age is the time since the onset of star-formation and r is the 
e-folding timescale for an exponentially decreasing SFH. is the stel- 
lar metallicity and E(B — V) is the colour excess which characterizes the 
Calzetti dust attenuation law. 



3 MODELLING HIERARCHICAL GALAXY 
FORMATION 

In this section we provide a brief description of the aspects of the 
GALFORM semi-analytic model which are relevant to this work. An 
introduction to the model and the associated underlying physics can 
be found m |Cole et al.] < [2000| i, |Baugh| p006| i and |Benson| ( [2l)IO] >. 
Briefly, GALFORM connects the properties of galaxies to a given 
cosmological model by coupling dark matter halo merger trees to 
a set of continuity equations that describe the exchange of baryons 
accreted onto dark matter haloes between stellar, cold disk gas and 
hot halo gas components. The physical processes that determine 
the form of these continuity equations include shock heating and 
subsequent radiative cooling of accreted gas onto galaxy disks, qui- 
escent star-formation in galaxy disks, chemical enrichment of the 
ISM, the ejection of cold gas and metals by supernova feedback, 
the suppression of gas cooling by AGN and photoionization feed- 
back and disk instabilities and galaxy mergers that can trigger both 
spheroid formation and bursts of star-formation. It is important to 
note that various versions of GALFORM have appeared in the liter- 
ature which we refer to as separate models. These models are dis- 
tinct from each other in that they all use different choices for model 
parameters and in some cases actually include physical processes 
which do not appear in other models. 



3.1 The Lagosl2 and Laceyl3 models 



We adopt the recently developed model described in Lag os et al.| 
( |2012} (hereafter Lagosl2) as the fiducial model to explore in this 
study. The Lagos 12 model is the most recent version of the model 
described in La gos et al.| (201 which in turn is a development 
of the model originally described in |Bower et al.|p006| > (hereafter 
Bower06). The Bower06 model was the first variant of GALFORM 
to include the effects of AGN feedback shutting down gas cool- 
ing in massive haloes. The Lagosl2 model is distinct from the 
Bower06 model in that it includes an alternative star formation 
law for galaxy disks based on an empirical relationship connecting 
the star formation rate in a galaxy to the molecular-phase gas den- 
sity. The molecular gas fraction is, in turn, related to the mid-plane 



pressure in the galaxy disk (Blitz & Rosolowsky 2006 1. This new 
law is observationally motivated and is characterized by parame- 
ters which are constrained by observations, greatly reducing the 
available parameter space in the model. Other changes with respect 
to Bower06 include longer timescales for both the minimum and 
total duration of starbursts and different reionization parameters. 
These changes were made to reconcile the model predictions with 
observations of LBGs (Lacey et al. 201 1 1. The model uses SPS files 
from a 1999 private release of the Bruzual & Chariot model family 
(BC99) and assumes a Kennicutt IMF ( Kennicutt 1983 ). Compared 
to the x = 1.35 Salpeter IMF used in the SED fitting procedure, 
the Kennicutt IMF has the same mass range (m* £ 0.1, 100 Mq) 
but has a steeper slope of x = 1.5 and a break in the power law 
at m* = IMq, below which the slope is x = 0.4. The IMF slope 
x is defined by d ^"^ = m~ x . The lack of a power law break 
at 0.1 Mq means that the Salpeter IMF has an overabundance of 
low mass stars compared to the Kennicutt IMF, resulting in higher 
M /L ratios. 

To help explore if certain aspects of our results are model 
dependent, we also consider the model presented in Lacey et al. 
(2013, in prep) (hereafter Laceyl3). The Laceyl3 model is a hy- 
brid of the Bower06 model family and the model from Bau gh et^T] 
(2005 ). It includes AGN feedback, starburst events triggered by 
disk instabilities and galaxy mergers, the star formation law de- 
scribed in Lagos et al. ( 201 lj and a non-universal IMF. We choose 
this model as a comparison to the Lagosl2 model because of sev- 
eral differences between the models which are relevant for SED 
fitting. The non-universal IMF used in the Lacey 13 model consists 
of a Kennicutt IMF for star formation in disks and a top-heavy IMF 
with slope x = 1 in starbursts. It should be noted that the slope 
of the top-heavy IMF used in the Lacey 13 model is less extreme 
than the x = top-heavy IMF required in Bau gh et al.| to match 
the number counts of submillimetre galaxies. The Laceyl3 model 
generates SEDs using the |Maras ton (2005) (MA05) SPS model. 
There has been a considerable amount of debate in the literature 
as to whether the luminosity of thermally-pulsating asymptotic gi- 
ant (TP-AGB) stars featured in the MA05 model is accurate (e.g. 
|Kriek et al.|2010| |Zibetti et al.|2013| >. This has potentially impor- 
tant consequences for the stellar mass inferred from SED fitting, 
potentially changing the M/L at NIR wavelengths by as much as 
50% for a stellar population of age w lGyr jMaraston et al.|2006| 
Michalowski et al. 2012). Although we do not explore this issue in 
any detail, the debate surrounding the contribution from TP-AGB 
stars makes the Lacey 13 model a useful comparison to our fiducial 
model. 



3.2 Calculating intrinsic galaxy SEDs 

The SED fitting procedure described in Section|2]relies on the accu- 
racy of SPS modelling to provide realistic SEDs for simple stellar 
populations. The same is true for GALFORM which uses SPS mod- 
elling to predict model galaxy SEDs from the star formation and 
chemical enrichment history of each galaxy, as calculated by the 
model. Compared to SED fitting, which has to assume a parametric 
form for galaxy SFHs and a single metallicity for all of the stars 
in a given galaxy, GALFORM self-consistently generates complex 
assembly histories for galaxies which include chemical evolution 
(see [Baugh 2006 for examples). There are only a small number 
of metallicities available for publicly released SPS models. There- 
fore, in order to actually use the chemical enrichment history of 
each galaxy in GALFORM, the model performs linear interpolation 
in \og(Z+) between the tabulated SSPs. This approach is not ap- 
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plied in standard SED fitting procedures which instead use a dis- 
crete metallicity grid. Another difference between GALFORM and 
SED fitting is that galaxies in GALFORM are divided into disk and 
bulge components. The net SED of each model galaxy predicted 
by GALFORM is therefore the sum of two composite stellar popu- 
lations, each with its separate star formation and chemical enrich- 
ment history. Finally, it should also be noted that GALFORM uses 
a different choice with respect to SED fitting regarding the treat- 
ment of the recycling of mass and metals. SPS models typically 
provide estimates of the amount of mass a SSP recycles back to the 
ISM as a function of age. This information is used in SED fitting to 
calculate the best-fitting stellar mass. For reasons of numerical effi- 
ciency, theoretical galaxy formation models, including GALFORM, 
typically do not use this information and instead apply the instanta- 
neous recycling approximation where mass and metals are instantly 
returned to the ISM. The amount of mass and metals returned per 
unit mass of stars formed are both fixed parameters in GALFORM 
and are therefore independent of the age of a galaxy. The exact re- 
cycled fractions and yields are calculated, for a given IMF, from the 
output of a SSP with solar metallicity at an age of 10 Gyr. This has 
a direct impact on how stellar mass is calculated in GALFORM and 
it is expected that this will lead to a small, redshift-dependent dis- 
agreement with the non-instantaneous recycling scheme employed 
in SED fitting. 

3.3 Dust attenuation 

The SED fitting procedure and GALFORM both use the same Madau 
( 1995 1 prescription for absorption and scattering of UV photons by 
neutral hydrogen in the IGM. However, there are important differ- 
ences in the way attenuation by dust is treated. The SED fitting 
procedure uses the Calzetti law which includes the assumption of 
a star-dust geometry corresponding to a uniform, foreground dust 
screen, as discussed in Section [272) GALFORM performs a more 
physical calculation of radiative transfer for a realistic geometry 
of the stars and dust. It models dust as a two phase medium sep- 
arated into diffuse dust in the ISM and compact dust clouds that 
enshroud star-forming regions ISilva et al.p998). For a d etailed 
description of this dust attenuation model, see |Lacey et al.|p011) 
and references therein and also Lacey, Baugh & Frenk (in prepa- 
ration). We provide a qualitative overview of the model here, fo- 
cusing on aspects of the modelling that are particularly relevant to 
our analysis. We use the standard terminology whereby extinction 
curves describe the absorption and scattering out of the sightline 
to a single star and attenuation curves describe the total absorption 
and scattering both into and out of all sightlines to an extended stel- 
lar distribution on the sky. We characterize attenuation curves with 
the effective optical depth, T e e, as a function of wavelength A. The 
effective optical depth, T e g,\, is defined by 

Toff, A = — ln(-F , atton,A/-Fintrin,A), (3) 

where -Fatten, a and i^ntrin.A are, respectively, the attenuated and 
intrinsic galaxy fluxes at a given wavelength and Fatten. a is calcu- 
lated from a radiative transfer model. 

The starting point for the dust model used in GALFORM is to 
calculate the total dust mass in each galaxy. This is calculated by 
assuming that the ratio of mass in dust to metals in the cold gas is 
a constant and that this ratio follows the value inferred for the lo- 
cal ISM (Savage & Mathis 1979). Dust is then divided into diffuse 
and compact birth cloud components. The relative fraction of dust 
mass in each component is a model parameter. The fraction in dif- 
fuse dust, /diffuse, is set to 0.75 in the Lagosl2 model and 0.5 in 



the Laceyl3 model. Both dust components use an input Milky-Way 
extinction curve as a starting point to calculate the resultant atten- 
uation curves of each galaxy by radiative transfer. It is important 
to realise that inclination and geometric effects can lead to total at- 
tenuation curves which are very different from the input extinction 
curve (e.g. Gonzalez- Perez et al.|2013) . 

The spatial distribution of diffuse dust is assumed to follow 
an exponential disk profile that traces the stellar disk, both in ra- 
dial and vertical scale-length. The effective optical depth associ- 
ated with diffuse dust is calculated by interpolating between the 
tabulated radiative transfer calculations performed by Ferra ra et al.| 
l |1999^ . [Ferrara et aT7| (p999 ) calculate the effective optical depth of 
disk-bulge systems as a function of wavelength, galaxy inclination, 
face-on extinction optical depth in the V-band, Tvo, and disk-to- 
bulge scale-length ratio, tvo is calculated directly from the density 
of dust at the centre of the disk, using the local ISM dust to metals 
ratio, and so scales with the surface density of diffuse dust in the 
galaxy disk as 

TVO OC /diffuse -McoldZcold/fdisk 2 , (4) 

where Afcotd and Z co \d are the mass and metallicity of cold gas in 
the galaxy disk and rdisk is the radius of the galaxy disk. It should 
be noted that no allowance is made in the modelling of diffuse dust 
for any differences in the relative spatial distributions of young and 
old stars in the galaxy disk, although this is accounted for with 
the second, compact dust cloud component. It is normally assumed 
that there is no dust in galaxy bulges. An exception is made, how- 
ever, for diffuse dust associated with gas forming stars in starbursts 
triggered by mergers or disk instabilities. In these starbursting sys- 
tems, the attenuation by diffuse dust is approximated by temporar- 
ily treating the bulge as a disk when using the results from the Fer- 
|rara et al~J (|T999 ) radiative transfer calculations. We discuss some of 
the advantages and potential problems associated with the way that 
diffuse dust is modelled in GALFORM, in the context of our results, 
in Section[6] 

The second dust component in GALFORM represents dust in 
dense molecular clouds enshrouding star-forming regions. As such, 
it generally affects the light emitted only by young stars which 
in turn are assumed to escape the dense dust clouds over a fixed 
timescale, which is a model parameter. The cloud component is 
therefore more significant in actively star-forming galaxies and 
starbursts where very young stellar populations can dominate large 
parts of the overall galaxy SED. The escape time is set to 1 Myr for 
the Lagosl2 model and 1 Myr for the Laceyl3 model. The clouds 
are modelled as being spherically symmetric with uniform density 
and a mass of 10 6 Mq and a radius of 16pc. The enshrouded stars 
are placed at the centre. Attenuation from this simple geometry can 
be evaluated analytically. For a more detailed description of this as- 
pect of the calculation, see Lacey, Baugh & Frenk (in preparation). 

The resultant combination of the diffuse and compact dust 
components attenuating the overall galaxy SED increases the level 
of physical realism beyond what is represented by the Calzetti law 
used in SED fitting. We pay particular attention to this in Sec- 
tion |4.5| but it should be noted that a full exploration of how the 
dust modelling used in GALFORM compares to the empirical rela- 
tions used in observational studies is beyond the scope of this paper. 
For more information, see Gonzalez- Perez et al.| (2013| l for a dis- 
cussion of how model predictions derived using this dust modelling 
approach compare with observations of LBGs. 
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4 STELLAR MASS RECOVERY 

In this section we first examine how accurately SED fitting can 
recover the stellar masses of a volume limited sample of model 
galaxies predicted by the Lagos 12 model at a selection of redshifts. 
We then attempt to isolate and explain the various different effects 
which affect the accuracy of stellar mass estimation. In this sec- 
tion, we use the standard SED fitting parameter grid and filter set 
described in the top half of Table [TJ 

4.1 Overview 

Fig. [TJ shows the ratio of the stellar mass estimated using SED 
fitting, M*[fit] to the true stellar mass in the Lagos 12 model, 
M+ [model], plotted as a function of M* [model] for a selection 
of redshifts. We choose to show individual galaxies colour coded 
by the density of points at a given position on the plane. We also 
show the 10, 50 and 90 percentile ranges of the distribution. This 
approach shows the broad trends in the overall distribution whilst 
still highlighting the presence of any unusual features or outliers. 
We quantify the distributions in each panel using two simple statis- 
tics in order to facilitate a rough quantitative comparison with other 
results presented in this section. We define /i as the mean value 
of the median offset in log 10 (M* [fit] /M* [model]) calculated for 
each bin in M* [model]. We define a as half of the mean value of 
the 68% range in Iog 10 (M* [fit] /M* [model]) calculated for each 
bin in M* [model]. If there is no dependence of the scatter and 
median offset on A/* [model], then /i and a quantify exactly the 
average systematic and random errors which affect the stellar mass 
estimation. 

At face value, the results shown in Fig. [Tj indicate that the ac- 
curacy of the stellar masses estimated using SED fitting is very 
poor, particularly at high redshift. It should be noted, however, 
that we have deliberately chosen to assume a Salpeter IMF in 
our SED fitting procedure despite the fact that the Lagos 12 model 
uses a Kennicutt IMF. The difference in M/L ratio between the 
Salpeter and Kennicutt IMFs can account for the systematic off- 
set in M* [fit] /Mi, [model] seen for low mass galaxies. However, 
the IMF mismatch cannot explain the behaviour displayed for mas- 
sive galaxies, particularly at high redshift. These galaxies display a 
huge scatter in Af* [fit]/M* [model] . Specifically, there seems to be 
a population of massive galaxies where the stellar mass is signifi- 
cantly underestimated. The medians and percentiles of the overall 
distribution show that this is an outlying population at low redshift. 
However, at high redshift, it is apparent that the stellar masses of 
almost all of the most massive model galaxies is significantly un- 
derestimated. In the most extreme individual cases, the stellar mass 
can be underestimated by factors greater than a hundred. Finally, 
the distributions display a level of bimodal behaviour which can be 
seen by eye from the point density distribution indicated by the 
colour scheme. This is easier to see in the higher redshift pan- 
els. The two peaks of the bimodal feature are typically offset in 
log 10 (Af* [fit] /Af* [model]) by « 0.25 dex. This is significant and 
clearly undesirable. 

There are a number of different factors in the SED fitting cal- 
culation that could combine to produce the behaviour shown in 
Fig. [TJ Therefore, it is useful to modify the SED fitting procedure 
and GALFORM in order to isolate how each factor of the calcula- 
tion contributes to this overall behaviour. In Fig. [2] we show how 
the distribution of estimated over true stellar mass changes with the 
inclusion or exclusion of these factors for the Lagosl2 model. We 
adopt a fiducial redshift of z = 2 for this exercise. Fig. [2^ shows 



z 0.5 1 2 3 4 

GALFORM: BC03 SPS, Salpeter / SED fitting: BC03 SPS, Salpeter, IRA 
fi/dex -0.01 -0.01 -0.03 -0.03 -0.01 -0.01 
cr/dcx 0.01 0.03 0.05 0.06 0.05 0.05 

GALFORM: BC03 SPS, Salpeter / SED fitting: BC03 SPS, Salpeter, NIRA 
fi/dex -0.01 -0.01 -0.01 0.01 0.04 0.05 
cr/dcx 0.01 0.03 0.05 0.05 0.05 0.05 

GALFORM: BC99 SPS, Salpeter / SED fitting: BC03 SPS, Salpeter, NIRA 
fi/dex -0.03 -0.01 -0.02 0.00 0.03 0.06 
cr/dcx 0.04 0.04 0.05 0.05 0.05 0.05 

GALFORM: BC99 SPS, Kennicutt / SED fitting: BC03 SPS, Salpeter, NIRA 
fi/dex 0.26 0.26 0.26 0.27 0.30 0.31 
cr/dcx 0.03 0.05 0.05 0.05 0.05 0.05 

Table 2. The mean median-offset fj, and half the mean 68% range, cr, of 
distributions in log(M* [fit]/Af* [model]) against M* [model] . All values 
listed are for the Lagos 12 model in the case where dust effects are ignored, 
both in the model and in the SED fitting procedure. Metallicity effects are 
also removed by forcing both the model and the SED fitting procedure to 
use Zi, = Zq. Each column corresponds to a different redshift. Each pair 
of rows corresponds to a different combination of choices made regarding 
SPS modelling, the IMF and recycling in GALFORM and the SED fitting 
procedure. The top pair of rows corresponds to the simplified case where 
GALFORM and the SED fitting procedure both use BC03 SPS models, in- 
stantaneous recycling (IRA) and a Salpeter IMF. The second pair of rows 
corresponds to the case where the SED fitting procedure is changed back to 
using default non-instantaneous recycling (NIRA). The third pair of rows 
corresponds to the case where the Lagos 12 model is changed back to using 
default BC99 SPS models and the SED fitting procedure uses NIRA. The 
final pair of rows corresponds to the default case where GALFORM uses 
BC99 SPS models, IRA and a Kennicutt IMF. The corresponding default 
SED fitting procedure uses BC03 SPS models, a Salpeter IMF and NIRA. 



the case where GALFORM and the SED fitting procedure have been 
stripped down to the point where effectively only the SFH is be- 
ing fit for each model galaxy. This is achieved by removing all of 
the effects associated with dust attenuation, chemical enrichment, 
recycling and the choice of SPS model and IMF. SPS model, re- 
cycling and IMF related effects are removed simply by making 
the two calculations consistent. Specifically, the SPS model used 
by GALFORM is changed to BC03 with a Salpeter IMF and the 
instantaneous recycling approximation is adopted in the SED fit- 
ting procedure. We remove chemical enrichment effects by forc- 
ing SED calculations in both GALFORM and the SED fitting pro- 
cedure to use solar metallicity. Dust effects are removed by setting 
E(B — V) = as a constraint in the SED fitting procedure and by 
using the unattenuated fluxes for model galaxies from GALFORM. 
From this simplified case, the other panels show how the distribu- 
tion changes with the reintroduction of the various aspects of the 
calculation that were removed in Fig. [2^. Each aspect is reintro- 
duced in isolation. 

The remainder of this section is outlined as follows. In Sec- 
tion [4^2] we discuss the role of SFHs, SPS models, recycling and 
the choice of IMF on the inferred stellar mass. In Section |4~3] we 
explore how our results are affected by wavelength coverage. Sec- 
tion |4.4| and Section |4"31 dis cuss the impact of metallicity and dust 
respectively. In Section |4T6] we extend our analysis to the Laceyl3 
model to explore the model dependence of our results. 
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Figure 1. The log of the ratio of the stellar mass estimated using SED fitting to the true stellar mass in the Lagosl2 model, plotted as a function of the true 
stellar mass. Each panel corresponds to a different redshift as labelled. The coloured points represent individual model galaxies. The point colours are scaled 
logarithmically with the local point density in the panel, from red at low density to yellow at high density. The black points and corresponding error bars 
show the median, 10 and 90 percentiles of the distribution in bins of true stellar mass, /i is the mean median-offset and a is half the mean 68% range of the 
distribution. For reference, the blue dashed line shows the locus of equality between estimated and true stellar mass. 



4.2 SFHs, recycling, SPS models and the IMF 

As discussed earlier, the case presented in Fig. [2^ is simplified to 
the extent where the only difference between SED calculations per- 
formed by GALFORM and the fitting procedure is in the form of the 



galaxy SFHs. The SED fitting procedure assumes an exponentially 
declining SFH characterized by the time since the onset of star- 
formation, i agc , and the e-folding timescale, r, whereas GALFORM 
self-consistently calculates the SFH of each model galaxy. None 
of the concerning features and trends seen in Fig. [T] are present in 
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Figure 2. The log of the ratio in the stellar mass estimated using SED fitting to the true stellar mass in the Lagos 12 model at z = 2, plotted as a function of 
the true stellar mass. Formatting of points and symbols is the same as in Fig.[T] The different panels show the distribution for different variations of both SED 
fitting and GALFORM. a) The SFHs of galaxies are the only factor which can vary in the SED fitting process. No dust extinction is applied to model galaxy 
SEDs in GALFORM and E(B — V) = is applied as a constraint in the SED fitting procedure. Z+ = Zq is applied as a constraint in both GALFORM and the 
SED fitting. The SPS model used by GALFORM is changed to BC03 with a Salpeter IMF (in order to be consistent with the SED fitting) and the instantaneous 
recycling approximation is used in the SED fitting procedure (to be consistent with GALFORM). b) SFHs, recycling, SPS models and the IMF are the only 
factors in the SED fitting process. Dust and metallicity related effects are removed as in Panel a), c) SFHs and metallicity are the only factors in the SED fitting 
process. Dust, recycling, SPS models and IMF related effects are removed as in Panel a), d) SFHs and dust are the only factors in the SED fitting process. 
Metallicity, recycling, SPS model and IMF related effects are removed as in Panel a). 



Fig. [2^, which instead shows a smooth distribution with a small 
scatter almost centered around the locus of equality between es- 
timated and true stellar mass. The distribution can be completely 
characterized by the mean offset /i = —0.03 dex and the mean 
spread a — 0.06 dex in this idealized case. It is perhaps surprising 
that, on average, the SED fitting works so well given the diversity 
of SFHs which can be predicted in GALFORM, and it is interest- 
ing then to see whether this result is reproduced at other redshifts. 
We list the fi and a values for this simplified case for other red- 
shifts in the top pair of rows in Table [2] It is interesting to see that, 
averaged over the entire galaxy population, the assumption of an 
exponentially declining SFH has almost no impact on the accu- 
racy of the stellar mass estimation at z — (fi — —0.01 dex and 
a — 0.01 dex). In addition, the small scatter seen in Fig. |2|i at 
z — 2 does not increase for higher redshifts. Comparing this level 
of scatter with that seen in Fig. [T] implies that, for our analysis, 
the assumption of an exponentially declining SFH has a negligible 
impact on stellar mass estimation. 



We strongly emphasize, however, that this result is only true 
for the average over all model galaxies with Af* 1.4 x 10 8 Mq. 
Closer inspection of Fig. [2^ reveals that there are outliers to the 
overall distribution where the stellar mass is underestimated by al- 
most an order of magnitude. In Fig. [3] we demonstrate that if a 
subset of the overall galaxy population is considered, the assump- 
tion of an exponentially declining SFH can lead to larger errors in 
the stellar mass estimation. In this case we compare the average 
offset and scatter of the entire galaxy population from the Lagos 12 
model at z — 2 with galaxies selected as being dominated by bursts 
of star-formation. We define bursts as galaxies with higher SFRs 
in a burst component relative to the SFR associated with quies- 
cent star-formation in galaxy disks. We choose this subset because 
these galaxies are likely to have SFHs that differ significantly in 
many cases from an exponentially declining SFH. Our definition 
of a burst is a straightforward physical definition that can be made 
in a theoretical model, and should not be confused with the typi- 
cal observational definition of a starburst as an object with an el- 
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Figure 3. The log of the ratio of the stellar mass estimated using SED fitting 
to the true stellar mass in the Lagos 1 2 model at z = 2, plotted as a function 
of the true stellar mass. As in Fig.[2ji, both GALFORM and the SED fitting 
procedure have been modified such that all dust and metallicity effects are 
removed and the IMF, SPS model and treatment of recycling are consistent 
between the two calculations. The top panel shows the distribution for all 
galaxies. The bottom panel shows the distribution for bursting galaxies, se- 
lected as model galaxies with a higher SFR in a burst component relative to 
the quiescent SFR in the galaxy disk. Formatting of points and symbols is 
the same as in Fig.[T| 



evated SFR compared to the mean of the population. Comparison 
of the distributions shown in Fig.|3]for all galaxies (top panel) and 
for bursts (bottom panel) shows that the stellar masses of bursting 
galaxies is underestimated on average by Afi — —0.1 dex com- 
pared to the average over the total galaxy population. In addition 
there is significantly increased scatter in M* [fit] / M* [model] when 
only bursts are considered. 

Fig. [2j> shows a similar scenario to Fig. [2^ but where the 
SPS model, IMF and treatment of recycling are changed back to 
be consistent with the default Lagos 12 model and default SED 
fitting procedure. Specifically, the Lagos 12 model uses BC99 
SPS models, instantaneous recycling and a Kennicutt IMF in this 
panel. The SED fitting procedure instead uses BC03 SPS mod- 
els, non-instantaneous recycling and a Salpeter IMF. Comparison 
of Fig. |2ji and [2|> shows that the different choices of SPS model 
and IMF, as well as the treatment of recycling, that can be made 
within SED fitting and GALFORM can cause constant offsets in 
[fit]/M* [model] but do not create additional scatter in the dis- 



tribution. It is interesting to explore the relative contribution from 
these different factors to these offsets. The remainder of Table [2] 
shows n and a over a set of redshifts for different combinations 
of choices regarding recycling, the SPS model and the IMF. In 
all cases, dust and metallicity effects are removed from both GAL- 
FORM and the SED fitting procedure. 

By default, the SED fitting procedure uses non-instantaneous 
recycling whereas GALFORM uses instantaneous recycling with a 
constant recycled fraction. This constant recycled fraction is fixed, 
for a given IMF, to the recycled fraction of a SSP with solar metal- 
licity and age 10 Gyr. There are two factors that could lead to sys- 
tematic errors in stellar mass estimation caused by differences be- 
tween instantaneous and non-instantaneous recycling. Firstly, the 
adopted relationship between initial and remnant mass for stars 
may be different in GALFORM and the BC03 SPS model used in 
the SED fitting procedure. [^] We check this by comparing the re- 
cycled fraction at 10 Gyr for a solar metallicity BC03 SSP with 
Salpeter IMF with the corresponding recycled fraction used by 
GALFORM for a Salpeter IMF. The two recycled fractions at 10 Gyr 
are R = 0.31 for the BC03 SSP and R = 0.30 for GALFORM 
which are almost consistent. We therefore do not expect this fac- 
tor to significantly affect the stellar mass estimation. Secondly, 
for non-instantaneous recycling, the recycled fraction is a func- 
tion of galaxy SFHs, whereas for instantaneous recycling the re- 
cycled fraction is independent of galaxy SFHs. Therefore, as the 
overall age distribution of the model galaxy population evolves, it 
is to be expected that part of any systematic error in stellar mass 
estimation caused by differences between instantaneous and non- 
instantaneous recycling will be redshift dependent. This is verified 
by comparing the values of /j, shown in the top and second sec- 
tions of Table [2] Changing from instantaneous recycling (top) to 
non-instantaneous recycling (second) in the SED fitting procedure 
has negligible impact at z — but results in a 15% offset in stellar 
mass by z = 4. This is a small effect compared to some of the other 
potential sources of error (e.g. the choice of IMF) but should still be 
accounted for if an attempt is made to make a precise comparison 
between stellar masses derived from observations and theoretical 
models that use instantaneous recycling, particularly at high red- 
shift. 

It is beyond the scope of this study to attempt to provide a 
comprehensive investigation into how stellar mass estimation is af- 
fected by uncertainties associated with SPS modelling and the IMF. 
Comparing the third and bottom sections of Table [2] shows that the 
difference between using a Salpeter and Kennicutt IMF is given by 
Afi ~ 0.25 — 0.29 dex for z 6 0,4. This simply demonstrates the 
well known result that changing from Salpeter to an IMF such as 
Kennicutt, that features a low-mass cutoff, results in M/L ratios 
that are offset by nearly a constant factor, reflecting the fact that 
stars at the low-mass end contribute a negligible amount to the in- 
tegrated light of a SSP. Comparing the second and third sections of 
Table|2]shows that there is also a small, ~ 7% increase in the scatter 
of the distribution at z — when the transition is made from us- 
ing BC03 SPS in GALFORM (second) back to the BC99 SPS model 
(third) used in the default version of the Lagosl2 model. This is 
presumably caused by some difference between these SPS models. 
However, given that the BC99 and BC03 SPS models both belong 
to the same overall model family, it should be noted that this is 



1 In GALFORM, we use the relations between initial and remnant masses 
from |Marigo"eraL]{T996) and |Portinari et al.|p998) . See |Cole et aTlpOOO) 
for details. 
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No NIR or IRAC filters - Observer Frame 
H/Aex -0.01 -0.02 -0.04 -0.06 -0.07 -0.07 
cr/dex 0.02 0.04 0.09 0.19 0.29 0.23 



Table 3. The mean median-offset fj, and half the mean 68% range cr of dis- 
tributions in M-i, [fit] I Mi, [model] against M*[model]. As in Fig. [2^, all 
values listed are for the Lagos 12 model in the idealized case where dust 
and metallicity effects are removed and the choice of SPS model, IMF and 
treatment of recycling is consistent between GALFORM and the SED fitting 
procedure. Each column corresponds to a different redshift. Each pair of 
rows corresponds to a different configuration of filters used to perform the 
SED fitting. The top pair of rows corresponds to the default case where the 
12 broad-band, observer-frame filters listed in Table [T] are used, spanning 
from -B435 to the 8.0/xm Spitzer IRAC band. The second pair of rows cor- 
responds to the same filter set with the modification that the filters are fixed 
in the galaxy rest-frame. The third pair of rows corresponds to a reduced set 
of observer-frame filters where the Spitzer IRAC filters are removed. The 
final pair of rows extends this by removing the J, H and K bands along 
with the IRAC filters. 



almost certainly a significant underestimate of the true impact on 
stellar mass estimation associated with uncertainties in SPS mod- 
elling. 



4.3 Wavelength coverage 

In the top section of Table [2] we show the mean offset /1 and mean 
spread a for a selection of redshifts in the idealized case where the 
Lagos 12 model and the SED fitting routine are stripped back to the 
point where only the SFH is different between the SED calcula- 
tions. For this idealized case, there is no systematic redshift depen- 
dence in the mean offset /i. However, there is a gradual increase in 
the scatter from a ~ 2% at z — up to a « 15% at z = 2. The 
scatter does not continue to increase beyond z — 2. The increase 
in scatter over the redshift range z £ 0,1 could be attributed to 
two separate effects. Firstly, any changes in the overall distribution 
of model galaxy SFHs with redshift could affect the accuracy of 
SED fitting. Over this redshift interval, there is a substantial fall 
with time in the overall SF activity, which may correspondingly re- 
flect a change in the underlying distribution of SFHs. Secondly, the 
rest-frame wavelength coverage of the filter set changes with red- 
shift such that the longer wavelengths in the rest-frame galaxy SED 
are no longer available in the SED fitting process at high redshift. 
To separate these two effects, we consider a modification of both 
the SED fitting procedure and GALFORM to use a filter set that is 
fixed in the galaxy rest-frame, independent of redshift. The top and 
second sections of Table [3] show values of fj, and mean spread a 
for the observer frame and rest-frame filter sets respectively. Using 
the rest-frame filter set removes most of the dependence of a on 
redshift, revealing that averaged over the entire galaxy population, 
any change in galaxy SFHs with redshift has a negligible impact on 



the accuracy of SED fitting when estimating stellar masses, at least 
when dust and chemical enrichment effects are not present. 

For the sake of completeness, it is also interesting to explore 
how important the NIR filters are for accurately estimating stel- 
lar mass in this idealized case where dust, metallicity, SPS model 
and IMF related effects have all been removed. The bottom four 
rows of Table [5] show fj, and a in the case where either the IRAC 
filters or all of the NIR filters are removed from the SED fitting 
process. Comparison to the full observer-frame filter set shown in 
the top part of Table [5] shows that having (perfect) photometry for 
the IRAC bands reduces the scatter in the stellar mass estimates by 
Act « 0.05 dex for z ^ 2. In the scenario where only the opti- 
cal filters are available, the accuracy of SED fitting degrades dra- 
matically above z = 1, even for the idealized scenario presented 
here. There is also an apparent trend whereby stellar masses are 
increasingly underestimated with increasing redshift. The degrada- 
tion at higher redshifts demonstrates that it is necessary to sample 
the rest-frame optical-NIR part of the intrinsic galaxy SED in or- 
der to properly account for the contribution from older stars which 
typically dominate the stellar masses of galaxies. 



4.4 Metallicity 

Fig.[2j: reintroduces metallicity variation back into GALFORM and 
the SED fitting procedure. For this panel, metallicity is a free pa- 
rameter in the SED fitting procedure and the full chemical enrich- 
ment histories of model galaxies are used to calculate their SEDs 
in GALFORM. Fig. [2J; demonstrates that the bimodal features seen 
in Fig. [T] are caused by some aspect of the SED fitting calculation 
associated with metallicity. To better understand this behaviour, we 
show in the top and bottom panels of Fig. [4] the same distribution 
with galaxies colour coded by their mean mass-weighted metallic- 
ity in GALFORM or by the best-fitting metallicity calculated in the 
SED fitting procedure. This reveals that while the metallicity of 
galaxies in GALFORM is continuous across the bimodal feature, the 
metallicity returned by the SED fitting procedure clearly traces the 
bimodality seen in Fig. [2};. This suggests that the SED fitting pro- 
cedure could be incorrectly associating a metallicity with a model 
galaxy because of degeneracies with other parameters such as age 
(e.g. |Bell & de J ong 20011. Underestimating the metallicity can 
lead to a corresponding overestimate of the age and hence the M/L 
ratio. Additionally, it is possible that the parameter grid of metal- 
licities used in SED fitting has insufficient resolution to reproduce 
the SEDs of model galaxies with mass-weighted metallicities that 
lie in between the values on the parameter grid. 

In order to understand what is causing the bimodal behaviour 
and to see if it can be removed, we have explored a number of 
different choices regarding how metallicity is treated in the SED 
fitting procedure. In Fig. [5] we show how these choices affect the 
distribution of estimated to true stellar mass against stellar mass 
for the Lagos 12 model at z — 2, for the case where dust, recycling, 
SPS model and IMF related effects have been removed. The first 
and most simple option we explore is simply to fix the metallicity 
of all galaxies to a constant value in the SED fitting procedure. 
This choice is often made in observational studies presented in 
the literature (e.g. |Rodighiero et al.|2010||Marchesini et al.|2009| >. 
The distribution for this case is shown in Fig. |5j). Although fixing 
the metallicity removes the bimodal behaviour, it also introduces a 
mass dependent bias into M* [fit] /M* [model], whereby the stellar 
masses of less massive galaxies is underestimated. This behaviour 
is clearly undesirable, although the problem might be alleviated 
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Figure 5. The log of the ratio of the stellar mass mass estimated using SED fitting to the true stellar mass in the Lagos 12 model at z = 2, plotted as a function 
of the true stellar mass. As in Fig.[2j;, both GALFORM and the SED fitting procedure have been modified such that all dust effects are removed and the IMF, 
SPS model and treatment of recycling are consistent between the two calculations. Formatting of points and symbols is the same as in Fig.^ Each panel 
corresponds to a different configuration of the SED fitting procedure, a) The default case as shown in Fig.|4]where metallicity is a free parameter in the fit and 
the mode (best-fit template) of the likelihood distribution is used to estimate the M / L ratio of each model galaxy, b) Metallicity is constrained to = Zq in 
the fit. c) Metallicity is forced to the closest possible value to the true mass-weighted metallicity of each model galaxy, d) Metallicity is a free parameter in the 
fit and the mean over the likelihood distribution is used to estimate the M/ L ratio of each model galaxy. The mean is calculated using a likelihood- weighted 
summation ^\ exp(— Xi 2 /2)M over the parameter space, e) Metallicity is a free parameter in the fit and additional template SEDs are added to the template 
grid by interpolating in metallicity. /) Metallicity is a free parameter in the fit, additional template SEDs are added to the template grid by interpolating in 
metallicity and the mean over the likelihood distribution is used to estimate the M/ L ratio of each model galaxy. 



somewhat if a restricted range in stellar mass is considered, as is 
often the case for high redshift galaxy samples. 

In Fig.|5J;, we force the SED fitting procedure to choose the 
closest metallicity on the parameter grid to the true mass-weighted 
metallicity of each model galaxy. This choice could only be repli- 
cated in an observational study if external constraints were avail- 
able on the stellar metallicity for each galaxy in the sample. Com- 
parison of the distribution shown in Fig. [3J: to the default case 
shown in Fig.[5^ shows that constraining the metallicity in this way 
restricts the bimodal behaviour to a narrow range in [model]. 
This in turn indicates that there is a degeneracy between two pos- 
sible metallicities which is broken when an external constraint is 
introduced. However, there is still a strong bimodality in the distri- 
bution at K'h [model] «4x 10 9 M Q . 

Another choice that can be made in SED fitting is to change 
the statistical method used to obtain the best estimate stellar mass. 
Instead of picking the point in the parameter space with the smallest 



X 2 (which corresponds to the mode of the likelihood distribution), 
it is also possible to estimate the stellar mass of each galaxy by cal- 
culating the mean over the likelihood distribution. This is achieved 
by performing a likelihood weighted summation £\ exp(— x* 2 /2) 
over the parameter space. Taylor et al. (2011) describe the imple- 
mentation and advantages of this weighted-average method in more 
detail. In principle, taking the mean rather than the mode should 
result in estimated stellar masses that are more robust against dis- 
creteness in the parameter space and could therefore help to re- 
move some of the bimodal behaviour seen in Fig. [5^ and Fig. |5j;. 
We show the distribution when this modified approach is used in 
Fig. BH, Comparison to Fig. [5^ reveals that the weighted-average 
approach does blur the bimodal feature, although the mass estima- 
tion is clearly still not perfect. 

As a final step, it is also possible to simply interpolate tem- 
plate SEDs between the metallicity points on the SPS metallicity 
grid. This will help especially in situations where the likelihood 
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Figure 4. The log of the ratio of the stellar mass estimated using SED fitting 
to the true stellar mass in the Lagos 1 2 model at z = 2, plotted as a function 
of the true stellar mass. As in Fig.[2j:, both GALFORM and the SED fitting 
procedure have been modified such that all dust effects are removed, and 
the IMF, SPS model and treatment of recycling are consistent between the 
two calculations. Top: Each point represents an individual galaxy from the 
Lagos 12 model and is coloured according the mean stellar mass- weighted 
metallicity calculated by GALFORM for that galaxy. Bottom: Each point rep- 
resents an individual galaxy from the Lagos 12 model and is coloured ac- 
cording to the best-fitting metallicity solution calculated in the SED fitting 
procedure. The parameter grid in available to the SED fitting procedure 
is shown in Table^ 



distribution associated with a single metallicity grid point is signif- 
icantly offset in M/L ratio from the others. In this case, neither the 
mean nor mode of the distribution will return a robust estimate of 
the true stellar mass if the outlying metallicity grid point dominates 
the overall distribution. We have confirmed that this is indeed the 
case for individual galaxies that fall on either side of the bimodal 
feature seen in Fig. [5^. We add points to the parameter space at 
Zi, — 0.5, 0.3, 0.1, 0.05 and 0.04 Zq using linear interpolation of 
the template SEDs in log(Z*). We show the effect of including this 
extended parameter grid in isolation in Fig.[5£ and combined with 
the weighted-average method in Fig. [5^. Adding the extra metal- 
licity points in isolation breaks up the main bimodal feature into 
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Figure 6. The log of the ratio of the stellar mass estimated using SED fitting 
to the true stellar mass in the Lagos 12 model at z = 2, plotted as a func- 
tion of the true stellar mass. As in Fig. |2ji, both GALFORM and the SED 
fitting procedure have been modified such that all metallicity effects are re- 
moved and the IMF, SPS model and treatment of recycling are consistent 
between the two calculations. Each point represents an individual galaxy 
and is coloured by the rest-frame V-band effective optical depth, Ty. c ff , as 
calculated in the Lagosl2 model. The colour scaling is indicated by the key. 

several smaller, less obvious features, improving the accuracy in 
the estimated stellar mass. Combining the interpolated metallic- 
ity grid with the weighted-average method almost completely re- 
moves any artificial bimodality in the distribution of the ratio of es- 
timated to true stellar mass against stellar mass. The overall scatter 
is also slightly reduced in the process. This approach could easily 
be adopted in observational SED fitting. We discuss this further in 
Section Wl2\ 

4.5 Dust attenuation 

In Fig.|2jl, we reintroduce dust attenuation back into GALFORM and 
the SED fitting procedure. Fig.|2ji shows that some aspect related to 
how dust attenuation changes galaxy SEDs results in a population 
of galaxies which are intrinsically massive but have stellar masses 
which are significantly underestimated by SED fitting. Fig.|6]shows 
the same distribution shown in Fig. |2ji but with individual galax- 
ies plotted as points coloured by their rest-frame effective optical 
depth in the V-band Tv, e ff , as calculated in the Lagosl2 model. This 
shows a clear trend whereby the SED fitting procedure systemati- 
cally underestimates the stellar masses of the model galaxies with 
the most dust extinction. Variations between the M/L ratio of tem- 
plate SEDs with sensible combinations of parameters are typically 
much smaller than the offset in M* [fit] /A/* [model] seen for these 
highly extincted galaxies. This indicates that rather than the prob- 
lem being caused by parameter degeneracies between e.g. dust and 
age, it seems that the SED fitting is simply not correctly recovering 
the overall normalization of the intrinsic model galaxy SED. This 
implies in turn that the Calzetti law must be a poor match to the net 
attenuation curves calculated for dusty galaxies in GALFORM. We 
confirm that this is indeed the case in Fig.|7]and Fig.[8] Fig.|7]shows 
the intrinsic and observed SEDs for 9 individual galaxies which are 
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Figure 8. Effective optical depth r e ff, plotted as a function of rest-frame 
wavelength at z = 2. Solid coloured lines show the attenuation of the La- 
gos^ model for the 6 galaxies shown in the middle and bottom rows of 
Fig. [7] The dashed grey lines show the attenuation curve from the Calzetti 
law for a range of values of E(B — V). The red dashed line corresponds to 
the specific case where E(B — V) = 1.0, which is the maximum reddening 
considered in the SED fitting procedure. 

selected to show a range of offsets in M* [fit]/M* [model] . It is im- 
mediately apparent that the SED fitting procedure underestimates 
the overall normalization of the intrinsic model galaxy SEDs for 
the dustiest galaxies. It is also particularly noticeable that the ra- 
diative transfer calculation used in the Lagos 12 model is applying 
a significant dust extinction to the entire SED, including up to the 
NIR for these galaxies. This behaviour cannot be reproduced by the 
Calzetti law. 

This can be seen more clearly in Fig. [8] which shows the atten- 
uation curves for the 6 model galaxies shown in the middle and bot- 
tom rows of Fig. [7] Also plotted is the Calzetti law for a wide range 
of values of E(B — V). The red dashed line corresponds to the 
Calzetti law with E(B — V) = 1, the maximum value we include 
in the parameter space used in the SED fitting procedure. Fig. [3] 
demonstrates that although the Calzetti law can match the amount 
of attenuation applied to the dustiest galaxies by the GALFORM ra- 
diative transfer model in the UV, it cannot reproduce the levels of 
attenuation at optical to NIR wavelengths. Furthermore, even if a 
different dust attenuation law was used in SED fitting that had the 
freedom to represent the types of attenuation curve of the model 
galaxies seen in Fig.[8] there would still be a very obvious degener- 
acy between the presence of a grey dust extinction component and 
simply having fewer stars producing light in a given galaxy. If such 
attenuation curves exist in reality, it would be very challenging to 
accurately constrain the stellar masses of dusty galaxies without 
performing detailed radiative transfer calculations using the entire 
UV - FIR SED. 

In order to help understand why some of the dusty galaxies 
in the Lagosl2 model have such extreme attenuation curves, it is 
also useful to consider the physical properties of the model galax- 
ies shown in Fig. [7] As discussed in Section[3] galaxies are divided 
in GALFORM between a disk and bulge component. In normal situa- 
tions, dust is assumed only to be present in the disk but this changes 
during starbursts, triggered by mergers and disk instabilities. Star- 
bursts are assumed to take place inside the bulges of galaxies and 
therefore a bulge dust component is required to account for dif- 
fuse dust in these systems. We find that the dusty galaxies where 



the SED fitting procedure fails are typically compact and have high 
star formation rates. The most extreme examples shown in the bot- 
tom row of Fig. [7] are compact starbursts. The intermediate cases 
shown in the middle row of Fig.[7]are compact, star-forming disks. 
It is beyond the scope of this paper to provide a detailed analysis 
of the relationship between the physical properties of model galax- 
ies from GALFORM and their dust attenuation curves. However, it 
is straightforward to see from Equation [4] that compact, gas rich 
galaxies will have the largest dust attenuations. In Section [63) we 
discuss the origin of the shape of the attenuation curves shown in 
Fig. [8] We also discuss whether observed galaxies with high dust 
content could have similar dust attenuation curves. 

4.6 Alternative galaxy formation models 

We explained in Section[T|that the focus for this study is not to at- 
tempt to provide an exhaustive, quantitative guide on the accuracy 
of stellar mass estimates derived from broad-band SED fitting. Part 
of the problem with using a SAM for this purpose is that results 
of this type will depend to some extent on the specific choices and 
assumptions made as part of that SAM. Instead, we focus on a spe- 
cific test case and try to understand and explain the origin of differ- 
ent systematics that appear in the distributions shown in Fig.[T] As 
an extension of this analysis, we also explore whether the behaviour 
seen in Fig. [T] is unique to the Lagos 12 model. In Fig. [9] we com- 
pare the distribution in estimated to true stellar mass against stellar 
mass from the Lagos 12 model at z — 2 with that in the Laceyl3 
model introduced in Section[3] We use the same SED fitting proce- 
dure used to fit model galaxies from the Lagos 12 model shown in 
Fig.fT] The two distributions shown in the top and bottom panels of 
Fig. Mare quite similar in some respects but notably different in oth- 
ers. The general trend whereby the stellar masses of progressively 
more massive galaxies is increasingly underestimated is seen for 
both models. However, unlike for the results for Lagos 12 model, 
the trend traced by the medians of the distribution continues mono- 
tonically all the way to the highest mass bins in the Lacey 13 model, 
implying that, on average, even the very most massive galaxies in 
the Lacey 13 model are very dusty. This implies in turn that they 
are forming stars at an elevated rate. This is an example of how 
our results can depend on the underlying physics that, in this case, 
controls the relative fraction of star-forming to passive galaxies. 

The a values calculated for the two distributions show that 
there is a slightly smaller level of random errors when estimating 
the stellar masses of galaxies from the Lagosl2 model, as compared 
to the Lacey 13 model. Visually inspecting percentiles reveals that 
this difference can be attributed to the smaller scatter in the dis- 
tribution in the low mass bins shown for the Lagosl2 model. The 
A/x ~ 0.3 dex offset between the two distributions can be under- 
stood as a product of two separate effects. Firstly, the recycled frac- 
tion associated with the Kennicutt IMF used in the Lagos 12 model 
is set to R = 0.39. The corresponding recycled fraction used for 
stars forming in disks in the Lacey 13 model is R — 0.44. This 
can account for A[i w 0.04 dex of the total systematic offset. Sec- 
ondly, the Lacey 13 model uses the |Maraston| i |2005| > (MA05) SPS 
model to compute galaxy SEDs whereas both the Lagos 12 model 
and the SED fitting procedure use versions of the Bruzual & Char- 
lot SPS model family. As discussed in Section [5] it is established 
that changing from using MA05 to BC03 SPS models in SED fit- 
ting of observed galaxies can change the estimated stellar masses 
by « 50 - 60% ( |Maraston et al.|2006||MTchalowski et al.|2012f . If 
the TP-AGB contribution is really as uncertain as the discrepancy 
between the BC03 and MA05 SPS models, the A/Lt » 0.3 dex off- 
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Figure 7. SEDs plotted as a function of wavelength in the observer frame for 9 model galaxies generated by the Lagosl2 model at z = 2. The black and 
blue solid lines show the dust attenuated and intrinsic SEDs respectively for the best-fitting SED templates calculated by the SED fitting procedure. The blue 
and red filled points show the intrinsic and attenuated flux respectively for model galaxies in each of the 12 photometric bands used in the fitting process. For 
plotting purposes, the SEDs are normalized such that the maximum flux of each model galaxy, as calculated by GALFORM, is zero. The three galaxy SEDs 
shown in the top row are selected quasi-randomly as examples where the SED fitting succeeds in recovering the intrinsic stellar mass. The three galaxy SEDs 
shown in the bottom row are selected as those with the biggest mass offsets in M* [Fit] / M t [Model] . The remaining three galaxy SEDs shown in the middle 
row are intermediate cases between these two extremes, log M* [Fit]/Af* [Model] is the log of the ratio of the estimated to true stellar mass. Ty. c ff [model] 
is the effective optical depth in the rest-frame V band, as calculated in the Lagosl2 model, Ty,cff [fit] is the effective optical depth in the rest-frame V band, 
as estimated by the SED fitting procedure. 



set between the two distributions has to be considered as a lower 
limit on the systematic uncertainty on stellar masses contributed by 
uncertainties from SPS modelling. 

Combined with the uncertainty in the IMF, these differences 
actually result in a total offset of /i = 0.53 for the Laceyl3 model, 
relative to the SED fitting using a Salpeter IMF. As discussed in 
Section[5] the Laceyl3 model uses a top-heavy IMF of slope x = 1 
in starbursts. We have checked whether this is important for stellar 
mass estimation by performing SED fitting on a modified version 
of the Laceyl3 model that uses a universal Kennicutt IMF. The dis- 
tribution for this scenario is shown in the middle panel of Fig. [9] We 
find that, averaged over the entire galaxy population, the distribu- 
tions with and without the top-heavy IMF are very similar, at least 
at 2 = 2. This is not unexpected because the optical-NIR SEDs 
(which to first order set the estimated stellar mass) of typical galax- 
ies are unlikely, on average, to be dominated by light directly emit- 
ted by starbursting populations. However, this will not necessarily 
be true for UV/FIR selected galaxy samples, particularly at higher 



redshifts. There is a small A/i « 0.04 dex mean offset between 
the two distributions, whereby the estimated stellar mass is slightly 
higher, on average, for the universal IMF version of the Laceyl3 
model. Visual inspection of the percentiles of the two distributions 
shows that the systematic associated with dust, seen at the high 
mass end in Fig. |2ji, is slightly less prominent for this modified 
version of the Lacey 13 model. The offset could therefore be caused 
by the higher rates of metal injection into the ISM that results from 
a top-heavy IMF in bursts. This in turn increases the dust content 
of the ISM in bursting systems, consequently increasing the impact 
of the systematics associated with dust. 

Finally, the bimodal behaviour seen in M* [fit] /M* [model] 
for the lower to intermediate mass bins in the Lagos 12 model is not 
apparent for the Lacey 13 model. We have investigated this further 
by performing SED fitting on the Lacey 13 model in the case where 
dust effects are ignored. In this case, we find that galaxies which 
are fit with different metallicities are, on average, systematically 
separated in recovered stellar mass. This is in agreement with the 
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Figure 9. The log of the ratio of the stellar mass estimated using SED fitting 
to the true stellar mass for different GALFORM models at z 2, plotted 
as a function of the true stellar mass. The top panel shows the distribu- 
tion for model galaxies from the default Lacey 13 model. The middle panel 
shows the distribution for model galaxies from a version of the Lacey 13 
model modified to use a universal Kennicutt IMF. As a reference, the bottom 
panel shows the distribution for model galaxies from the Lagos 12 model, as 
shown in Fig.^ Formatting of points and symbols is the same as in Fig.^ 
In all cases, the SED fitting uses a Salpeter IMF 



results seen in Fig.Elfor the Lagosl2 model. However, this effect 
is not visible in Fig7j9] because any underlying bimodality in the 
distribution is blurred out by the larger overall scatter at low masses 
seen for the Lacey 13 model. 



5 THE STELLAR MASS FUNCTION 

In Section [4] we show that there are various systematics that can 
prevent the SED fitting procedure, described in Section[2] from ac- 
curately estimating the stellar masses of individual model galaxies 
calculated by different GALFORM models. We now turn our atten- 
tion to addressing the question of how these systematics affect the 
global statistics of the galaxy population. We explore this issue by 
comparing the intrinsic stellar mass functions predicted by different 
GALFORM models to the corresponding mass functions recovered 
from the same models using SED fitting. Measurements of the stel- 
lar mass function are often used to constrain hierarchical galaxy 
formation models (e.g. |Guo et al.|20fT ; Henrique s et al.|2012) . It is 
therefore useful to also make a comparison with different observa- 
tional estimates of the stellar mass function. This also helps to put 
the systematic effects explored in Section|4]into context. 

We present stellar mass functions for a selection of redshifts 
from the Lagosl2 model in Fig.[l0]and from the Laceyl3 model in 
Fig. [TT| We also present stellar mass functions for the GALFORM 
model described in Baugh et al. ( 2005 ) in an appendix. We modify 
our standard SED fitting configuration (outlined in Table [TJ at this 
point by assuming a Chabrier IMF instead of a Salpeter IMF. This 
choice is made in order to be consistent with the bulk of the obser- 
vational studies shown in Fig. [To] and Fig.[TT] For simplicity, we 
do not choose to change our SED fitting procedure between each 
redshift panel. It should be noted, however, that the different obser- 
vational studies all use slightly different variations of SED fitting 
parameters and filter sets. In addition, the low redshift observational 
studies ( |Li & White|2009[|BaId~ry et al.|2012) use alternative SED 
fitting methods compared to the standard procedure described in 
Section|2] Baldry et al. ( 2012 1 use the likelihood-weighted summa- 
tion technique, as explored briefly in Section |4~4| and described in 
detail in |Taylor et al.| < |20 1 1 ^ . [Li & WhiteH2009) use stellar masses 
calculated with the nonnegative matrix factorization technique de- 
scribed in Blanton & Roweis (2007). For consistency with the re- 
sults shown at other redshifts, we do not use these alternative meth- 
ods for the mass functions recovered from the model. 

Comparison of the intrinsic model mass function (solid-blue 
line) with the mass function recovered using SED fitting (solid-red 
line) in Fig. |10| shows that the systematics seen in Fig. [TJcan have 
an appreciable impact on the inferred global statistical properties 
of the galaxy population. The intrinsic and recovered model mass 
functions agree best at the low mass end but disagree at the knee of 
the mass function and at the high mass end. This becomes increas- 
ingly evident in the higher redshift panels. The dominant factor re- 
sponsible for this disagreement is dust, as discussed in Section pO] 
This is demonstrated by comparing the recovered stellar mass func- 
tions when dust effects are (solid-red line) and are not (dashed-red 
line) included in the Lagos 12 model and the SED fitting process. 
The recovered stellar mass function that includes dust extinction 
effectively cuts away the knee of the intrinsic model mass function 
for z > 0. This is not seen when dust effects are not included, 
which can be understood by comparing Fig. [2^ with Fig. [2jl. A 
common feature of the recovered mass functions, both including 
and not including dust, is that in the highest redshift panels, the 
abundance of the most massive galaxies is increased with respect to 
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Figure 10. Stellar mass functions predicted by the Lagos 12 model for a selection of redshifts, as labelled in each panel. The solid blue line shows the intrinsic 
stellar mass function produced by the Lagos 12 model. The solid red line shows the stellar mass function recovered using SED fitting when dust effects are 
included and a Chabrier IMF is assumed in the fitting procedure. As a reference, the dashed red line shows the corresponding stellar mass function where no 
dust extinction is applied to the model galaxy SEDs and E(B — V) = is used as a constraint in the fitting procedure. The grey points and error bars show 
observational estimates of the stellar mass function from |Li & Whi te ( 2009 1, Baldry et al. 1 2012 1, llbert et al. ( 2010"), |Santini et al.|j2012) and |Mortlock et al.| 
(201 1) . Where necessary we convert these observational results from a Salpeter to a Chabrier IMF using a —0.24 dex correction, calculated by comparing the 
recovered stellar mass using Salpeter and Chabrier IMFs with BC03 SPS models. 
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Figure 11. Stellar mass functions predicted by the Lacey 13 model for a selection of redshifts, as labelled in each panel. The definition and formatting of the 
lines and points is the same as in Fig. |10| 



the intrinsic stellar mass function predicted by the model. This sim- 
ply reflects the Eddington bias where the exponential decline of the 
mass function at the massive end means that any scatter in the stel- 
lar mass estimation shifts more galaxies into higher mass bins than 
vice- versa. This effect competes with the impact of dust attenuation 
to give the resultant shape of the solid red line in Fig.[T0] Compari- 
son to the observational data shows that, in some cases, the relative 
differences between the intrinsic and recovered model stellar mass 



functions are much smaller than the disagreement with the obser- 
vational estimates of the stellar mass function, particularly at the 
low-mass end and at low redshift. In such cases, the model clearly 
does not accurately reproduce the observational data and the stel- 
lar mass function can be used as a meaningful constraint. However, 
in the higher redshift panels the uncertainties on the observational 
data are larger and the differences between the intrinsic and recov- 
ered model stellar mass functions also become larger. Taking the 



© 2013 RAS, MNRAS 000,[Tp7l 



The stellar masses of galaxies 19 



difference between the recovered and intrinsic model mass func- 
tions as a measure of the level of uncertainty regarding what model 
curve should actually be compared to the data, it is then apparent 
that it becomes difficult to place any meaningful constraints on the 
model using the observed stellar mass function at high redshift. 

Fig. [TT] shows the same information as Fig. [To] but for the 
Laceyl3 model instead of the Lagos 12 model. The relationship be- 
tween the intrinsic and recovered stellar mass function is similar to 
what is seen in Fig.[lO]but there are also a number of differences. 
The most obvious difference is that the overall density normaliza- 
tion of the recovered stellar mass function is higher than for the 
intrinsic stellar mass function. This occurs predominantly because 
of differences between the MA05 SPS model used in the Laceyl3 
model compared to the BC03 SPS model used in the SED fitting 
procedure. As discussed in Section |4~6| this means that the stellar 
masses of individual model galaxies are systematically overesti- 
mated for the Laceyl3 model by SED fitting, shifting the overall 
stellar mass function to the right. Comparison of the two recov- 
ered model stellar mass functions when dust effects are (solid red 
line) and are not (dashed red line) included shows behaviour sim- 
ilar to what is seen in Fig. [10] whereby the perceived abundance 
of massive galaxies at and above the knee of the mass function is 
suppressed, with the level of suppression increasing toward higher 
redshifts. It is particularly striking that the entire shape of the mass 
function can change dramatically, even in the lower to intermediate 
redshift panels. At these redshifts, the intrinsic stellar mass func- 
tion predicted by the Laceyl3 model does not resemble a single 
Schechter function as there is an apparent change in the power-law 
slope before the break. This feature is washed out in the recovered 
stellar mass function, which instead resembles a single Schechter 
function when dust effects are included. Finally, as an aside, it is 
interesting to note that there is a fairly strong level of disagreement 
between the shape of all of the model stellar mass functions and the 
shape of the observed mass function at z = 0. A disagreement in 
the overall shape cannot be explained by any systematic uncertainty 
such as the M/L associated with the IMF, and is interesting given 
that the model is tuned to reproduce the z = K-band luminosity 
function. This suggests that the stellar mass function can, at least at 
low redshift, provide useful constraints for galaxy-formation mod- 
els that are complementary to those provided by luminosity func- 
tion data. 

5.1 Lyman-break galaxies 

Up until this point we have focused on the relationship between the 
recovered and intrinsic stellar mass predicted for galaxies at low to 
intermediate redshift. At these redshifts, typically the photometric 
errors are small and colour selections do not need to be used in 
order to obtain well defined galaxy samples. Given that we find that 
the systematics in the distribution of log(M* [fit] /M* [model]) can 
have an appreciable impact on the recovered stellar mass function at 
these redshifts, it is interesting to explore how these effects translate 
to high redshift Lyman break galaxy (LBG) samples, where the 
photometric and redshift errors can become large, colour selections 
exclude parts of the galaxy population and the available optical to 
NIR wavelength coverage in the rest frame is reduced. 

We modify our SED fitting set-up at this point in order to more 
closely resemble the choices made by Stark et al. 12009|> and [Lee] 
|et al.| p012), who estimate the stellar mass function of LBGs. This 
alternative parameter grid is outlined in the bottom half of Table [T] 
and is described in Section [2] We also consider a number of ef- 
fects relevant for LBG samples that were not considered previously. 
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Figure 12. Stellar mass functions predicted by the Lagos 12 model at z = 4. 
The top panel demonstrates how the intrinsic stellar mass function predicted 
by the model is reshaped as a result of LBG selection criteria. The dashed 
blue line shows the intrinsic stellar mass function of all galaxies predicted 
by the Lagosl2 model. The dash-dotted blue line shows the intrinsic stellar 
mass function of galaxies selected using the B-dropout criterion from Stark 
|eTaTlp009) if no flux errors are included. The solid blue lines, present 
in both panels, show the corresponding intrinsic stellar mass functions of 
LBGs when the selection includes the effects of artificial flux errors. The 
model galaxy fluxes are artificially perturbed to mimic the S/N for each 
band quoted in Table 1 of Lee et al. 1 20121. The bottom panel demonstrates 
how the stellar mass function of LBGs, as recovered using SED fitting, is 
reshaped as a result of errors on the photometry and the statistical method 
used to construct the mass function. The red line shows the stellar mass 
function of LBGs recovered by SED fitting. The green line shows the cor- 
responding recovered mass function when the model galaxy fluxes are arti- 
ficially perturbed. This line lies underneath the red and magenta lines at all 
but low masses. The magenta line shows the corresponding recovered mass 
function when the model galaxy fluxes are artificially perturbed and the full 
stellar mass PDF of each individual galaxy is used to construct the mass 
function. The points and corresponding error bars show the observationally 
inferred stellar mass functions from|Stark et al. ( 2009) and |Lee et 31^2012) . 

Firstly, we explore the effect of artificially perturbing the fluxes of 
model galaxies, following a Gaussian distribution consistent with 
the 5cr limiting depths listed in Table 1 of |Lee et al.H2012} . These 
limiting depths are also used to calculate <r„ in Eq.[T]when perform- 
ing SED fitting. When the flux of a model galaxy falls below the la- 
limiting depth in a given band, we use the method of dealing with 
non-detections described in Section 3.1 of |Lee et al.| ( |2012) . Sec- 
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ondly, we explore the effect of LBG selection, taking the dropout 



criteria from Stark et al. 1 2009 1 (both S/N and colour selection cri- 



teria), which extends in redshift beyond [Lee et al.|j2012) to include 
z ~ 6 1775 dropouts. When exploring how the SED fitting performs 
without including the effects of dust attenuation, we still use atten- 
uated fluxes for the purposes of deciding which galaxies pass the 
LBG selection criteria and to decide which bands are counted as 
detections for individual galaxies. This ensures the galaxy samples 
are consistent when comparing different recovered mass functions. 
Finally, [Lee et al.| ( |2012} construct their observed stellar mass func- 
tions by summing together the stellar mass probability distribution 
functions (PDFs) calculated by their SED fitting procedure for each 
individual galaxy, hereafter referred to as the PDF method. This is 
distinct from standard practice in SED fitting where a single stel- 
lar mass value is assigned to each galaxy by considering only the 
best-fitting template to a given galaxy. This difference could poten- 
tially become significant at the low mass end of the mass function 
where galaxies are sufficiently faint that they are only detected in 
filters sampling the rest-frame UV. UV photometry does not pro- 
vide strong constraints on the stellar mass associated with older 
stars in galaxies and this uncertainty will be accounted for on an 
object-by-object basis when using the PDF method to construct 
the stellar mass function. To provide a fair comparison to the mass 
functions from |Lee et al.|pO"T2"] >, we also explore the effect of using 
the PDF method on our recovered mass functions. 

In Fig.[T2] we demonstrate how these various choices and op- 
tional modifications affect the stellar mass function predicted by 
the Lagos 12 model at z — 4. The top panel shows how the intrin- 
sic mass function predicted by the Lagosl2 model is reshaped by 
LBG selection and artificial flux errors. LBG selection criteria are 
designed to isolate UV-bright, star-forming galaxies over a given 
redshift range. For the Lagos 12 model, this has the effect of re- 
ducing the normalization of the mass function in a fairly uniform 
manner over the range of masses considered here. The impact of in- 
cluding artificial flux errors has minimal impact on LBG selection. 
The bottom panel shows how the stellar mass function of LBGs, 
as recovered by SED fitting, is affected by artificial flux errors and 
the PDF method. All of the recovered mass functions shown are 
virtually identical apart from at the low mass end. This simply re- 
flects the fact that more massive galaxies are typically brighter and 
are therefore detected with higher overall S/N and at optical-NIR 
wavelengths where the photometry is not as deep. 

In Fig. [T3] we present recovered and intrinsic model stellar 
mass functions from the Lagos 12 and Laceyl3 models for a se- 
lection of redshifts. All of the model mass functions shown are 
constructed using LBG selection criteria and artificially perturbed 
model galaxy fluxes. The PDF method is used to construct the re- 
covered mass functions, consistent with |Lee et al.|pOT2) . We also 
show measurements of the mass functions of LBG-selected sam- 

and |Lee et al.| (2012) . It should be 



Stark et al. 



Stark et al. 



2009 



(2009) apply an absolute magnitude cut at 



pies from 
noted that 

A4i500 = ^20, the effect of which can be clearly seen as the abun- 
dance of galaxies starts to fall below Af* w 2 x 1O 9 M0 . Stark et al. 
(2009) also give their results for a Salpeter IMF which we correct 
by —0.24 dex. This correction was estimated by comparing the re- 
covered stellar mass using Salpeter and Chabrier IMFs with BC03 
SPS models. When comparing the model mass functions with ob- 
servational data, it should be noted that we do not attempt to mimic 
errors associated with photometric redshifts for our model galax- 
ies. Including redshift errors is likely to increase the abundance of 
galaxies at the high-mass end, in line with observational samples 



that do not attempt to account for the Eddington bias associated 
with redshift errors. 

On first examination, Fig. [13] shows a similar picture to that 
seen in Fig. [To] and Fig. [TT] The inclusion of dust effects in the 
SED fitting procedure reshapes the recovered stellar mass function 
at the intermediate to high mass end (as seen by comparing the solid 
and dashed red lines). It is interesting to see that while the intrinsic 
model stellar mass functions (blue lines) are quite a good match to 
the observed stellar mass functions in many cases, the model mass 
functions recovered using SED fitting are in very poor agreement 
when dust effects are included. This emphasizes the danger of di- 
rectly comparing intrinsic mass functions predicted by theoretical 
models to observational data at high redshift, without accounting 
for the relevant uncertainties. It is also notable that the significant 
differences between recovered mass functions that include (solid 
red lines) and do not include dust effects (dashed red lines) are 
present in all of the panels shown. This demonstrates that for the 
models considered in our analysis, dust continues to play a role in 
reshaping the entire UV-NIR SEDs of massive galaxies all the way 
out to z = 6. Overall, similar to the situation seen in Fig. |10|and 
Fig. [TJJ it is striking that the relationship between the recovered 
and intrinsic model stellar mass functions shown in Fig. [T3] can 
vary dramatically, dependent on the different SED fitting choices 
which are made and the redshift considered. 



6 DISCUSSION 

6.1 The role of the SFH in estimating accurate stellar masses 

Estimating the physical properties of galaxies from SED fitting re- 
quires the adoption of a prior distribution of galaxy SFHs. Large 
variations exist in the prior distributions used in different obser- 
vational studies, reflecting the overall uncertainty in the optimal 
choice of SFH distribution. As a consequence, considerable effort 
has gone into establishing how these choices can affect the different 
galaxy properties that can be estimated from SED fitting (e.g. [Lee 
let al.|2009UMaraston et al.|2010||Pfbrr et al.|201^|Michatowski 
|et al.||2012| |Banerji et al.|[2013] [Schaerer et al.|[20T3] ). |Lee et al. 
(2009) and Schaerer et al. (2013) arrive at similar conclusions in 



that they both find that the stellar masses of LBGs are not particu- 
larly sensitive to the choice of SFH prior. At face value, this seems 
to conflict with the conclusions from other studies which find that 
the estimated stellar masses of specific galaxy classes at high red- 
shift are strongly sensitive to the assumed SFH distribution (e.g. 



|Maraston et al.|2010| [MichaIows ki et al.|2012||Banerji et al.|2013| >. 



It should be noted, however, that |Michalowski et al. (2012 ) and 
|Banerji et al.| ( |20i"3"] > allow for the possibility of multi-component 
SFHs. This choice gives the SED fitting procedure more freedom 
to fit a very young, UV-bright stellar population at the same time 
as including a significant population of older stars that contribute 
to the total SED primarily at longer wavelengths. This approach is 
shown to yield systematically higher stellar masses relative to as- 
suming a smooth exponentially declining SFH or an instantaneous 
starburst. 

We have demonstrated in the top section of Table [2] that, av- 
eraged over the entire galaxy population at a given redshift, the 
assumption of an exponentially declining SFH does not lead to a 
mean systematic offset in stellar mass. This is consistent with the 
results from |Lee et al.| ( [2009| ) who also find that, averaged over a 
population of model LBGs, stellar masses are well constrained us- 
ing an exponentially declining SFH. In addition, we have shown 
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Figure 13. Stellar mass functions predicted by the Lagosl2 (left side) and Laceyl3 (right side) models for a selection of redshifts, as labelled in each panel. 
The dashed blue lines show the intrinsic mass functions predicted by the models without imposing any selection criteria. In all of the other cases, model galaxy 
samples are constructed using the LBG selection criteria from |Stark et al.|(2009| . The flux of each galaxy is artificially perturbed to mimic the S/N for each 
band quoted in Table 1 of Lee et al. (2012). The solid blue lines show the intrinsic stellar mass functions of LBG-selected galaxies predicted by GALFORM. 
The solid red lines show the recovered model stellar mass functions of LBG-selected galaxies when dust effects are included in the SED fitting procedure. The 
dashed red lines show the corresponding recovered model stellar mass functions when dust effects are removed after selection criteria have been applied. All 
of the recovered mass functions are constructed using the full stellar mass PDF of each individual galaxy. The points and corresponding error bars show the 
stellar mass functions from Stark et al. (2009) and |Lee et aflpOLiY 
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that the scatter associated with fitting the SFH of model galaxies 
is negligible relative to the random and systematic errors caused 
by other factors in the SED modelling. It should be noted that this 
result only applies strictly to the case of fitting the SFHs of model 
galaxies in isolation, independent of dust and chemical enrichment 
effects. In practice, assuming a given prior SFH distribution may 
play a larger role in creating errors in stellar mass estimates because 
of the degeneracies that exist between the effects of age, metallic - 
ity and dust. For example, Pforr et aL] ( [2012| > show that in some 
cases it can be advantageous to ignore dust entirely when apply- 
ing SED fitting to galaxies with SEDs dominated by older stellar 
populations. This prevents the SED fitting procedure from incor- 
rectly fitting young and highly reddened galaxy templates to these 
galaxies. 

We have also demonstrated in Fig. [5] that if galaxies are se- 
lected as being burst dominated, the errors on the estimated stel- 
lar masses associated with assuming an exponentially declining 
SFH are greatly increased relative to the average error for the total 
galaxy population. This is not a particularly surprising result; we 
have simply selected galaxies for which a single-component, expo- 
nentially declining SFH is least likely to be appropriate. However, it 
does serve to reconcile our results with the findings of Michalowski 
et al. (2012) in the sense that fitting smoothly varying SFHs to 
galaxies that arc expected to have undergone recent bursts of star- 
formation (e.g. submillimeter galaxies) can lead to systematically 
underestimating the stellar masses of these objects. 

6.2 How should metallicity be included in SED fitting? 

It is standard practice to either fix the metallicity or to use a 
small number of discretely spaced metallicities when performing 
broad-band SED fitting to estimate the physical properties of galax- 
ies. Typically, interpolation between the metallicities is not used. 
This situation can be attributed to a combination of the inabil- 
ity of broad-band SED fitting to constrain metallicity (e.g. |Pforr| 
|et al.|20 12), the small number of SSP metallicities made available 
for popular SPS models and for reasons of numerical efficiency. 
Our analysis has shown that using a discrete and sparsely sam- 
pled metallicity grid causes undesirable bimodal features in the 
distribution of log(M*[fit]/M*[model]), as seen over a specific 
range in stellar mass in Fig. [4] We have also shown in Fig. 
that the decision to instead fix the metallicity of all galaxies to 
= Zq leads to an equally undesirable mass-dependent bias in 
log(M* [fit] I Mi, [model]). This behaviour is straightforward to re- 
move. Fig. [5J demonstrates that using interpolation to add more 
metallicities to the parameter grid as well as taking the mean rather 
than the mode of the probability distribution to estimate stellar mass 
can resolve any biases in stellar mass estimation associated with 
metallicity. This works primarily because interpolation acts to fill 
the gaps between metallicity grid points that are significantly offset 
in M/L ratio for a given age but also because because taking the 
mean helps to blur out discreteness in the distribution of M/L ra- 
tios for a given parameter grid (as discussed in |Taylor et al.|201 1| >. 
These steps could easily be incorporated into the standard SED fit- 
ting procedures used in observational studies. 

The success of the revised fitting procedure shown in Fig. [5j 
suggests that the standard assumption of a single stellar metallicity 
for all of the stars in a galaxy is acceptable, provided that interpo- 
lation and the method of performing a likelihood-weighted aver- 
age over all templates is used. This is perhaps surprising given that 
model galaxies in GALFORM can have complex chemical enrich- 
ment histories. In particular, the metallicity of the individual stellar 



populations that make up model galaxies is strongly correlated with 
age. |Conroy et al.H2009| > find that there is no significant difference 
between the optical and NIR colours of a single metallicity SSP and 
those of a multi-metallicity SSP of the same average metallicity. 
This would imply that SED fitting should correctly infer the M/L 
ratios of galaxies, provided that the procedure selects the correct av- 
erage metallicity for a given galaxy. This is entirely consistent with 
our results. However, |Conroy et al . ( 2009 1 note that they do not ac- 
count for a correlation between metallicity and age in their analysis 
which could feasibly create a significant difference in the colours 
of single and multi-metallicity stellar populations. Given that such 
correlations exist for model galaxies in GALFORM, it is therefore 
reassuring that Fig.]^ shows this effect does not have a significant 
impact on the estimated stellar masses of galaxies. |Gallazzi & Bell| 
( 2009 1 do find that fitting single or double colours of mock galaxies 
with chemical enrichment histories that contain an age-metallicity 
correlation has an impact on inferred M/L ratios. However, they 
find that this effect is small, consistent with our results. It should be 
noted that there is still a small level of mass-dependent bias evident 
in Fig.|5f. Given the results of |Gallazzi & Bell|p009"l >, this could 
potentially be explained as a result of fitting model galaxies' multi- 
metallicity stellar populations with single metallicity templates. 

Stellar mass and stellar metallicity are extremely strongly cor- 
related in the Lagosl2 model, as can be seen in Fig. [4] This is 
potentially significant in the context of our analysis because the 
strength of this correlation may serve to exaggerate the strength 
of the bimodal features seen in Fig. [4] In addition, if the stellar 
metallicity of real galaxies at a fixed stellar mass differs from that 
of model galaxies predicted by the Lagos 12 model, then the mass 
scale where any bimodal behaviour appears will change, compared 
to the feature seen in Fig. [4] We also note that it is specifically the 
lowest sub-solar metallicities which are responsible for the offsets 
seen in Fig. [4] If the stellar metallicities of real galaxies, at a given 
stellar mass, are higher than predicted by the Lagos 12 model, then 
it is possible that these sub-solar metallicities will not be relevant 
for the galaxies probed in most observational samples. In this case, 
the size of the discreteness effects seen in Fig. [4] will be signifi- 
cantly reduced. 

6.3 Can galaxies have significant dust attenuation at optical 
to NIR wavelengths? 

The most significant source of error we encounter in estimating the 
stellar masses of model galaxies is found when SED fitting is ap- 
plied to very dusty model galaxies. In Fig. [8] we have shown that 
these galaxies have much larger amounts of attenuation at longer 
wavelengths than is possible from the Calzetti law. This causes the 
stellar masses of dusty model galaxies to be significantly underes- 
timated. In extreme cases, this underestimate can be by factors as 
large as several hundred. From Fig. [I] it can be seen that this error 
affects the majority of model galaxies above 10 10 Mq at z > 2. 
Consequently, for these redshifts, the errors associated with dust 
can completely reshape the recovered stellar mass functions shown 
in Fig. [TO] Accordingly, the intrinsic rest-frame optical and NIR lu- 
minosity functions predicted by GALFORM will also be reshaped by 
dust. This will have a significant impact on any attempt to compare 
this particular model with observational data at higher redshifts. 

In order to understand this overall result, it is useful to con- 
sider how there can be significant dust attenuation at optical to NIR 
wavelengths for model galaxies in GALFORM. Firstly, it is impor- 
tant to note that for the assumption of a star-dust geometry corre- 
sponding to a uniform foreground dust screen, and for a specific 
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dust grain model, the ratio of absolute extinction Ay relative to 
the reddening E(B - V) must be constant. [Calzetti et al7| j2000) 
assumed this star-dust geometry and applied energy balance ar- 
guments to UV and FIR observations of 4 local starbursts to fix 
Rv = A V /E(B -V) = 4.0 5. Once fly has bee n fixed in this 
way, the attenuation curve from |Calzetti"etaLl (1994) 

can no longer 

reproduce the attenuation curves shown in Fig. [8] 

In reality, the assumption of a uniform dust screen is a very 
poor approximation to a realistic star-dust geometry. In GALFORM, 
disk stars are embedded in a diffuse dust component with the same 
spatial distribution as the stars. In this case, the path length through 
the diffuse dust to the observer will be different for each star. There- 
fore, the light from each star will experience a different amount of 
attenuation, yielding a net attenuation curve for the entire galaxy 
which can be significantly different from the input extinction curve 
of the dust grains (Gonzalez-Perez et al. 2013). A simple exam- 
ple that demonstrates this behaviour is to consider a star-dust ge- 
ometry corresponding to an infinite uniform slab containing stars 
and dust mixed together with the same uniform spatial distribution. 
The effective optical depth, r e a, for this geometrical configuration 
is given by 



T c ff,A 



hi 



1 — exp (— t ,a seci) 



r , a sec i 



(5) 



where i is the inclination angle of the slab relative to the observer 
and to,a is the face-on extinction optical depth for a single sight- 
line through the slab. In the limit that to,a becomes large, so that 
the slab is optically thick, this simplifies to T c ff,A — In (to, a seci). 
In this scenario, nearly all of the light emitted by stellar popula- 
tions from within the slab is absorbed and only light from a layer 
of stars at the surface can reach the observer. The light that escapes 
the slab only passes through a small amount of diffuse dust and 
as a consequence, is only reddened by a small amount. Therefore, 
for the optically thick case, this configuration will yield net atten- 
uation curves which, compared to the Calzetti law, are consider- 
ably greyer. This explains how the attenuation curves of very dusty 
galaxies in GALFORM can have significant amounts of attenuation 
at long wavelengths. 

In reality, dust in galaxies is not thought to exactly trace the 
spatial distribution of stars in galaxy disks. By comparing the red- 
dening of nebula emission (i.e. through the H a to Hp line luminos- 
ity ratio) with the reddening observed in the total UV continuum of 
galaxies, it is apparent that dust in the ISM must be concentrated 
around star-forming regions (e.g. Calzetti et al. 1994). This obser- 
vation has motivated the use of two-component dust models (e.g. 
|Silva et al.| 1 998 ) that contain a compact, birth-cloud dust compo- 
nent that attenuates the light emitted by very young stellar popu- 
lations. This scheme is applied in GALFORM and as such, the net 
attenuation curves of model galaxies shown in this paper take this 
effect into account. However, the diffuse dust in galaxies is also 
not thought to exactly trace the spatial distribution of stars. Specif- 
ically, the scale height of diffuse dust is known to be smaller than 
the overall scale height of stars in galaxy disks, particularly when 
considering older stellar populations (e.g. |Wild et al.|20TT) . 

Additionally, the presence of a clumpy ISM means that there 
are likely to be some sightlines through galaxy disks which are rel- 
atively free of dust. This effect was considered by Conroy et al. 
(2010b) who explored the effect of a clumpy ISM with a lognormal 
column density distribution of dust combined with the empirical 
dust model from |Charlot & Fall| < f2000| ). Although they show that 
the dumpiness of the ISM (characterized by the width of the log- 
normal column density distribution) can have a large impact upon 



NUV and optical colours, they state that there is negligible im- 
pact from dust on the 7f-band luminosity function. It should be 
noted that this conclusion depends strongly on how they calcu- 
late the total amount of dust in the diffuse ISM. If there is a suf- 
ficiently large mass in diffuse dust, a clumpy ISM will produce 
similar emergent behaviour in terms of the shape of the total at- 
tenuation curves to that of the slab geometry considered earlier. 
That is, the total SED of a dusty galaxy with a clumpy ISM will be 
dominated by light emitted by stars that lie along relatively unob- 
scured sightlines (which experience only minimal reddening) while 
light emitted from behind or within optically thick regions will be 
completely absorbed in comparison. However, the effect of having 
a clumpy ISM will, for an equal mass in diffuse dust, result in a 
lower normalization of the total attenuation curve as compared to 
an optically thick, uniform slab. This is simply because a higher 
fraction of the stars will be unobscured for the case of a clumpy 
ISM. At present, the dust modelling in GALFORM does not account 
for any dumpiness of diffuse dust in the ISM. This may result in a 
higher normalization for the attenuation curves of very dusty model 
galaxies, compared to real dusty galaxies. Although the inclusion 
of a birth-cloud dust component does account in part for dumpi- 
ness in the ISM, the impact from birth clouds will be of secondary 
importance if the diffuse dust component contains enough mass to 
absorb all of the light not emitted from close to the surface of the 
galaxy disk. 

The problem of choosing an appropriate star-dust geometry 
becomes even more complex in the case of a galaxy merger. It is 
well established from numerical simulations that pressure forces 
experienced by gas in the ISM can decouple the spatial distribu- 
tions of gas relative to stars, such that gas is funnelled into a com- 
pact region in the centre of the system, producing a nuclear burst 
of star-formation. Wuyts et al. (2009) show that by applying SED 
fitting to a suite of idealized hydrodynamical simulations of galaxy 
mergers, the stellar masses of simulated galaxies can be system- 
atically underestimated. This occurs because as the diffuse dust is 
concentrated into the central region, any light emitted from within 
or behind that region will be almost entirely cut out of the ob- 
served galaxy SED. In addition, the overall stellar distribution will 
be much more spatially extended than the gas during this phase and 
consequently will suffer minimal reddening. This is exactly analo- 
gous to the behaviour discussed previously for an embedded star- 
dust geometry or a clumpy ISM. However, as noted by Wuyts et al. 
(2009), the situation is complicated in this case by the fact that the 
stellar populations which are heavily obscured are younger, on av- 
erage, compared to the total stellar population. A strong correlation 
between stellar population age and the dust column density will 
serve to dilute the greying effect discussed for the slab and clumpy 
ISM examples (see Section 4.3.2 in |Wuyts et al.|2009| >. At present, 
GALFORM fails to account for these effects in specific situations; in 
the event of a major merger or disk instability, stars of all ages are 
mixed evenly with diffuse dust and gas in the galaxy bulge. This 
geometry is unlikely to be representative of real merging systems 
(although the situation is far less clear at high redshift) and as such, 
the attenuation curves of systems in GALFORM that are undergoing 
major mergers or disk instabilities may be unrealistic. In this case, 
the impact of dust on stellar mass estimation could be exaggerated 
to some extent for these systems. Changing the radial scale length 
of the burst compared to the stellar bulge would represent only a 
small change to the current implementation in GALFORM and we 
plan to investigate the impact of this change in future work. 

Aside from uncertainties associated with the star-dust geom- 
etry, it is important to appreciate the uncertainties associated with 
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calculating the mass and density of dust in the ISM of galaxies, 
particularly at high redshift. Generally speaking, theoretical galaxy 
formation models that attempt to model dust use local relations that 
give the ratio of dust to metals in the ISM ( |Cole et al.|20 00). These 
local relations are then applied universally, which is a large extrap- 
olation in the case of actively star-forming galaxies at high redshift 
where the physical conditions in the ISM can be very different. 
In addition, other aspects of a given theoretical galaxy formation 
model will have a strong impact on the final effect of dust on model 
galaxy SEDs. For example, in the case of our analysis, if the calcu- 
lations of metallicities or galaxy sizes are incorrect in GALFORM, 
then the size of the errors associated with dust on stellar mass esti- 
mation will also be incorrect. 

In summary, the assumption of a uniform foreground dust 
screen must be incorrect in most cases for real galaxies, but the 
details of the star-dust geometry in real galaxies may also be more 
complex than what is assumed in GALFORM. In addition, both the 
mass and density of dust calculated by GALFORM is dependent on 
both the overall accuracy of the model and being able to extrap- 
olate the local dust to metal ratio up to high redshift. We plan to 
explore how these factors could affect our results for stellar mass 
estimation in future work. 

Having explained why dusty galaxies in GALFORM can have 
significant amounts of attenuation at longer wavelengths relative to 
the Calzetti law, it is useful to consider if there is any evidence for 
this behaviour from other observational or theoretical studies. IPforrl 
|et al.|p0"T2| > and |Lee et al. ( 2009 1 do not find any evidence that dust 
can cause the stellar masses of galaxies to be significantly underes- 
timated when they fit model galaxies from other SAMs. However, 
we emphasize that the SAMs considered in these studies do not 
attempt to model dust attenuation in a physically motivated way. 
Instead, they adopt the same Calzetti attenuation curve in the SAM 
as in the SED fitting. It is therefore of no surprise that these authors 
do not recover the same results shown by our analysis. 

|Lo Faro et allpOBl l fit the full UV-FIR SEDs of an observed 
sample of 3 1 luminous and ultraluminous infrared galaxies at z = 1 
and 2 = 2. From their Figures 2, 5 and 7 it can be seen that their 
fitting method estimates significant attenuation of the stellar contin- 
uum across the entire UV-NIR SED for a number of objects in their 
sample. They also apply a standard UV-NIR SED fitting procedure 
to the same galaxy sample and find that, for the most dust obscured 
galaxies, the stellar mass obtained in this case is underestimated rel- 
ative to what is estimated from the full UV-FIR fitting procedure. In 
the top-left panel of their Figure 6, they show a clear trend whereby 
stellar mass is increasingly underestimated by UV-NIR SED fitting 
for increasingly dust obscured systems. This is in qualitative agree- 
ment with our results. It should be noted that they use the radiative 
transfer code GRASIL ( |Silva et al.|1998| l to generate UV-FIR tem- 
plates. GRASIL assumes the same star-dust geometry as is assumed 
in GALFORM. On the one hand, this assumed geometry is physi- 
cally motivated and is clearly superior to the crude assumption of a 
uniform foreground dust screen. However, the various uncertainties 
discussed earlier associated with choosing this particular star-dust 
geometry are also relevant for their analysis. 

|Michalowski et al.| ( |2010| > also used GRASIL to fit the full 
UV-FIR SEDs of a sample of 76 spectroscopically confirmed sub- 
millimeter galaxies, estimating stellar masses for these objects. 
Michalows ki et al.| {2012| > then revisited the same sample but in- 
stead applied a standard UV-NIR SED fitting procedure to estimate 
stellar masses. As discussed in Section [6~T| the intention behind 
their analysis was to investigate how priors on the SFH distribu- 
tion of submillimeter galaxies can affect stellar mass estimation 



for this class of objects. However, in the context of this discus- 
sion, it is interesting to consider the top-right panel of Figure 2 in 
Michalowski et al. ( 20121, where the stellar masses estimated using 
standard UV-NIR SED fitting are compared to the stellar masses 
calculated using GRASIL modelling of the full UV-NIR SED from 
|Michalowski et aLl pOlO). Submillimeter galaxies correspond to 
the objects in our analysis where the stellar mass would be most 
affected by optical-NIR attenuation. Therefore, it is striking that 
in contrast to the results from |Lo Faro et aT] ([20131, there does 
not appear to be a significant systematic difference between the 
stellar masses estimated using the standard UV-NIR and GRASIL- 
based UV-FIR methods of fitting submillimeter galaxy SEDs, at 
least compared to the uncertainties associated with choosing an ap- 
propriate SFH. 

It is important to note that the SED fitting method applied in 
|MichaIowski et al.|p0T0l > differs fro m|Lo Faro et al.|j2 013 1 in that 
they use set of template SEDs from[Tglesias-Paramo et al.|p007] >. 
These templates were constructed for a limited range of the pos- 
sible parameter space in GRASIL, chosen to reproduce the SEDs 
of star-forming galaxies in the local Universe. In contrast, Lo Faro 
|et al.| ( |2"013| > do not impose any strong priors on the various free 
parameters in GRASIL, exploring a large parameter space. This dif- 
ference in approach could potentially explain how stellar mass es- 
timation from standard UV-NIR SED fitting is only found to be 
strongly affected by |Lo Faro et al.| ( |2013) . Clearly, for the com- 
plex UV-FTR SED fitting procedures applied byfMichalowsk Tet al.| 
(2010) and |Lo Faro et al.] ( |2013} , the resultant stellar masses will 
depend strongly on the priors and assumptions that are adopted. 

|da Cunha et aLlpOTO] ) fit the UV-FIR SEDs of 16 local ultralu- 
minous infrared galaxies using an alternative procedure to GRASIL. 
Although they do not attempt to compare the stellar masses of these 
objects estimated using their UV-FIR SED fitting method with what 
would be estimated from standard UV-NIR SED fitting, it can be 
clearly seen that their fitting procedure favours a significant level of 
dust attenuation at optical-NIR wavelengths for all of the objects in 
their sample. This is consistent with the behaviour revealed by our 
analysis and with |Lo Faro et al.]^2013) . 

In the local Universe, the problems associated with assuming 
a specific star-dust geometry can be lessened for resolved, low in- 
clination galaxies. |Zibetti et al. (2009) show that when optical-NIR 
SED fitting is applied locally to derive the stellar mass density at 
each pixel in the images of resolved galaxies with prominent dust 
lanes, the total stellar mass calculated can be higher by up to 40 % 
relative to the stellar mass obtained by fitting integrated photome- 
try. This result is consistent with the trends shown in our analysis 
and can be understood as the result of the light emitted by stars that 
reside either within or behind optically thick dust lanes being sub- 
dominant in the total galaxy SED, compared to the light emitted by 
stars on unobscured sightlines. Unfortunately, this method cannot 
be readily extended to high redshift where very dusty galaxies are 
more common. 

6.4 How should theoretical galaxy formation models be 
compared to observational data? 

Our analysis is intended to demonstrate that aside from the well 
documented uncertainties on stellar mass estimation associated 
with SPS modelling and the form of the IMF (e.g. Conr oy et al.| 
2009 ), stellar mass estimation can also be significantly affected by 
the combined effects of dust, metallicity and recycling. However, 
viewed from another perspective our method of applying SED fit- 
ting to model galaxy SEDs offers, in principle, a new way to com- 
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pare predictions that involve stellar mass to the results from ob- 
servational studies. This approach is attractive because it allows 
models to be self-consistently compared with different observa- 
tional datasets without the need to change various parameters in the 
model (SPS model, IMF) in each instance to make a fair compar- 
ison. In addition, the scatter between intrinsic and estimated stel- 
lar mass is self-consistently accounted for using this method. This 
alleviates the need to invoke an arbitrary level of scatter in order 
to create agreement between model predictions and observations 
(e.g. |Guo et al.|20TT) |Bower et al.|2012) . In the case where pre- 
dictions from a theoretical model are compared simultaneously to 
both observables and inferred quantities such as stellar mass, it is 
clear that our methodology should be followed to make the com- 
parisons self-consistent with each other. Otherwise, the process of 
transforming from observables to intrinsic galaxy properties will 
become confused. For example, assumptions made in theoretical 
models to predict luminosity functions are likely to be in conflict 
with the assumptions made in SED fitting to estimate stellar mass 
functions. Unless our methodology is followed, using both of these 
diagnostics at the same time could therefore adversely affect any 
attempt to constrain the underlying physics of galaxy formation. 



It is important to realise that our methodology does not avoid 
the problem of converting intrinsic galaxy properties into observ- 
ables. Instead, the burden of accounting for these uncertainties is 
simply shifted from the observational SED fitting procedure back to 
the theoretical modelling process. The natural alternative is to only 
consider intrinsic galaxy properties and leave all of the uncertain- 
ties in the observational process of estimating these quantities. This 
approach is widely used in both the semi-analytic modelling and 



hydrodynamical simulation communities (e.g. |Bower et al. 


2006 


Neistein & Weinmann||2010| |Dave et al.||2011| |Guo et al. 


2011 


Khochfar & Silk||201 1| |Lagos et al.||201 1| |Lamastra et al. 


2013 


Ciambur et al. 2013]).[Conroy et al. ( 2010b I advocate using this lat- 



ter approach with the caveat that derived quantities such as stel- 
lar mass should only be calculated with the inclusion of a full 
marginalisation over all of the relevant uncertainties in SPS and 
dust modelling. Such a process would require that these uncertain- 
ties can be fully characterized by a discrete set of parameters and 
that robust prior distributions for the plausible ranges of these pa- 
rameters can be found. Our study serves to emphasize that choosing 
an appropriate distribution of priors would be extremely difficult. 
For example, the extinction optical depth of diffuse dust, assumed 
to be small in |Conroy et al.|pbl0b} , turns out to play a very impor- 
tant role in stellar mass estimation for our analysis when the optical 
depth becomes large. 

Finally, Fig. [10] and Fig. [TT] show that the systematics and 
biases that affect stellar mass estimation in our analysis cannot 
fully account for the level of disagreement between the stellar mass 
functions predicted by GALFORM models and the mass functions 
estimated from observational data. Given that these models were 
tuned to reproduce luminosity function data, this demonstrates that 
stellar mass functions contain complementary information to lumi- 
nosity functions, independent of the uncertainties associated with 
converting intrinsic galaxy properties into observables. This is en- 
couraging and suggests that provided the comparison between the- 
oretical models and observations is made self-consistently, existing 
estimates of the stellar mass function can provide significant con- 
straining power for galaxy formation models. 



7 SUMMARY 

Motivated by the desire to understand whether stellar mass is an ap- 
propriate tool for constraining hierarchical galaxy formation mod- 
els, we have used the observational technique of SED fitting to es- 
timate the stellar masses of model galaxies from the semi-analytic 
model GALFORM. Following the standard SED fitting procedure 
for fitting broad-band photometry, we find that effects associated 
with metallicity, recycling and dust can bias stellar mass estimates. 
In some specific cases, these effects can create systematic errors in 
stellar mass that are comparable to or greater than the potential sys- 
tematic errors associated with the uncertain form of the IMF. Fur- 
thermore, we have shown that these error sources are often stellar 
mass dependent, such that the stellar mass function of model galax- 
ies recovered using SED fitting can differ substantially in shape as 
well as in normalization from the intrinsic mass function predicted 
by a given model. 

The cause and nature of the individual systematic error 
sources uncovered by our analysis are as follows: 

• The exponentially declining star-formation histories that are 
typically assumed in SED fitting do not, averaged over the entire 
galaxy population, create any significant systematic errors in stellar 
mass. In addition, when averaged over the entire galaxy population, 
the random errors in stellar mass caused by fitting with exponen- 
tially declining SFHs are small. This is demonstrated in Fig. [2^ and 
Table [2] These results are selection dependent. If model galaxies 
are deliberately selected to be undergoing bursts of star formation, 
the assumption of an exponentially declining SFH leads to both 
systematic underestimates and a significantly larger scatter in the 
estimated stellar masses of these systems. 

• Differing assumptions regarding recycling of mass from stars 
back into the ISM can lead to small, redshift dependent system- 
atics in stellar mass. These are outlined in Table [2] Theoretical 
galaxy formation models typically apply the instantaneous recy- 
cling approximation, whereas standard SED fitting procedures use 
the best-fitting template SFH to estimate the recycled mass. Neither 
approach will give the correct answer in detail. Furthermore, the 
systematic differences between the two approaches should be ac- 
counted for when the stellar masses predicted by theoretical models 
that assume instantaneous recycling are compared to observational 
data. 

• Metallicity has the effect of introducing undesirable bimodal 
features into the distribution of recovered stellar mass that can be 
seen in Fig|2j;. This behaviour arises because the standard SED 
fitting procedure uses discrete, poorly sampled metallicity grids 
and a statistical method of choosing only a single best-fitting tem- 
plate (the mode of the probability distribution). Alternatively, if 
the equally common choice of fixing metallicity in SED fitting is 
implemented, the resultant estimated stellar mass suffers a strong 
mass-dependent bias. These problems can be solved in a straight- 
forward manner by following two simple steps. Firstly, interpola- 
tion can be used to fill in the gaps of the original metallicity grid 
provided for publicly available SPS models. Secondly, the statisti- 
cal technique advocated by Taylor et al. ( 201 1) can be implemented 
where the mean over the probability distribution, calculated from a 
likelihood-weighted average over all templates, is used to calculate 
a best estimate for the stellar mass of a given galaxy. 

• Dust attenuation in massive, dusty galaxies causes standard 
SED fitting procedures that assume a Calzetti law to systematically 
underestimate stellar mass. This occurs because the radiative trans- 
fer calculations performed in GALFORM predict significant attenua- 
tion at optical-NIR wavelengths in some cases. Thus, the light emit- 
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ted by obscured stellar populations in these galaxies is not properly 
accounted for when estimating stellar mass. 

We find that the shape of the stellar mass function at z — 
is robust against these error sources. However, at higher redshifts 
the systematic errors associated with dust significantly reshape the 
recovered mass functions such that a clear break in the intrinsic 
model mass function at these redshifts can be blurred out. We are 
forced to conclude that any attempt to constrain theoretical galaxy 
formation models using stellar mass functions from high redshift 
galaxy samples should only be performed with great care, given 
the potential for large mass-dependent systematics in stellar mass 
estimation from SED fitting. 
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APPENDIX: STELLAR MASS FUNCTIONS FOR THE 
BAUGH05 MODEL 

Throughout the main text we consider the Lagos 12 model from |La-| 
|gos et ak] ( |2012[ > and the Laceyl3 model from Lacey et al. (2013, 
in prep). In this appendix, we also consider the model described 
in |Baugh et al.| l[2005) (hereafter Baugh05). The Baugh05 model 
is distinct from the Lagos 12 and Lacey 13 models in that it does 
not include bursts of star formation triggered by disk instabilities 
or the updated star-formation law described in Lagos et al.| ( (20TT| ). 
The Baugh05 model also uses different timescales for star forma- 
tion, both in galaxy disks and in bursts, and uses supernova driven 
superwinds instead of AGN feedback as a mechanism to suppress 
the bright end of the luminosity function. Finally, as described in 
Section [3~T] the Baugh05 model uses a top-heavy IMF in bursts 
with a slope of x — 0. This is more extreme than the x = 1 slope 
used in the Lacey 13 model. 

We present stellar mass functions for a selection of redshifts 
from the Baugh05 model in Fig. [14] Neither the intrinsic or recov- 
ered stellar mass functions agree with the observational estimates 
of the stellar mass function at z = 0. The model overpredicts the 
abundance of low mass galaxies at z ^ 1 and overpredicts the 
abundance of the most massive galaxies at z — 0, suggesting that 
the feedback schemes implemented in this model could be unreal- 
istic. Similar behaviour regarding the effect of dust on the recov- 
ered stellar mass functions is seen with respect to the Lagos 12 and 
Lacey 13 models. At z = 0, the recovered stellar mass functions 
(both including and excluding dust attenuation effects) are lower in 
normalization with respect to the intrinsic model mass function. 
This could be a result of the SPS models used in the Baugh05 
model. Alternatively, the difference could be caused by the top- 
heavy IMF in bursts. We will explore this in more detail in future 
work. 
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Figure 14. Stellar mass functions predicted by the Baugh()5 model for a selection of redshifts, as labelled in each panel. The solid blue line shows the intrinsic 
stellar mass function produced by the Baugh05 model. The solid red line shows the stellar mass function recovered using SED fitting when dust effects are 
included and a Chabrier IMF is assumed in the fitting procedure. As a reference, the dashed red line shows the corresponding stellar mass function where no 
dust extinction is applied to the model galaxy SEDs and E(B — V) = is used as a constraint in the fitting procedure. The grey points and error bars show 
observational estimates of the stellar mass function from |Li & Whi te (20091, Baldry et al. 1 2012), [TTbert et al.|pQ10) , |Santini et al.H2012) and |Mortlock et al.| 
(201 1) . Where necessary we convert these observational results from a Salpeter to a Chabrier IMF using a —0.24 dex correction, calculated by comparing the 
recovered stellar mass using Salpeter and Chabrier IMFs with BC03 SPS models. 
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